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ABSTRACT 

Hydrogen Lya is our primary emission-line window into high redshift galaxies. 
Surprisingly, despite an extensive literature, Lya radiative transfer in the most realis- 
tic case of a dusty, multi-phase medium has not received detailed theoretical attention. 
We investigate Lya resonant scattering through an ensemble of dusty, moving, opti- 
cally thick gas clumps. We treat each clump as a scattering particle and use Monte 
Carlo simulations of surface scattering to quantify continuum and Lya surface scat- 
tering angles, absorption probabilities, and frequency redistribution, as a function of 
the gas dust content. This atomistic approach speeds up the simulations by many 
orders of magnitude, making possible calculations which are otherwise intractable. 
Our fitting formulae can be readily adapted for fast radiative transfer in numerical 
simulations. With these surface scattering results, we develop an analytic framework 
for estimating escape fractions and line widths as a function of gas geometry, mo- 
tion, and dust content. Our simple analytic model shows good agreement with full 
Monte-Carlo simulations. We show that the key geometric parameter is the average 
number of surface scatters for escape in the absence of absorption, Ao, and we provide 
fitting formulae for several geometries of astrophysical interest. We consider two inter- 
esting applications: (i) Equivalent widths. Lya can preferentially escape from a dusty 
multi-phase ISM if most of the dust lies in cold neutral clouds, which Lya photons 
cannot penetrate. This might explain the anomalously high EWs sometimes seen in 
high-redshift/submm sources, (ii) Multi-phase galactic outflows. We show the charac- 
teristic profile is asymmetric with a broad red tail, and relate the profile features to 
the outflow speed and gas geometry. Many future applications are envisaged. 



1 INTRODUCTION 



The hydrogen Lya line is our primary emission line window 
on the high-redshift universe. It is almost invariably cru- 
cial in securing red shift-identification s for the high e st red - 
shift galaxies ( e.g. .iHu et al.l j2002aiP: lAiiki et al.l J2003t); 
iKodaira et aP (l2003^ : iRhoads et alJ ll2003D : ISantos et all 
j200^V r. Besides yielding redshifts, the shape of the line 
profile, equivalent width, and offset from other emis- 
sion/absorption lines also encode information about the ge- 
ometry, kinematics and underlying stellar population of the 
host galaxy. For instance, features in Ly a emission have bee n 
used to suggest strong galactic outflows (Kunth ct al.'lQO^, 
as a signature of strong accretion shocks (Barkana fc Loeb, 
l2003l) . and as evidence for an unusuall y strong ionizing con- 
tinuu m, perhaps due to Pop III stars jMalhotra fc Rhoadj 
|2002|) . Even after escaping the environs of the host galaxy, 
Lya photons undergo processing in the surrounding inter- 
galactic medium (IGM), and the presence or absence of ob- 
served Lya emission can be used to place constraints on 
the epoch of reionization (iFan et al 1 12002) : iHaimanl 120021 : 
ISan tos 2004). Because of the numerous factors which con- 



tribute to Lya radiative transfer, the interpretation of such 
features is fraught with complexity. For instance, in a com- 
prehensive set of ~ 1000 Lyman Break galaxies at z ~ 3, a 
plethora of Lya strengths and line shapes were seen, ranging 
from pure damped absorp tion, to emission plu s absorption, 
to pure strong emission dShaplev et al 1 12003^ . Because of 
the tremendous potential returns for interpreting some rich 
data-sets, it is crucial to strive for a more detailed theoretical 
understanding of Lya emission line features. 

A very important factor in Lya transmission is the 
presence of dust. Since the massive stars which produce 
metals evolve on a short timescale, an d indeed superso- 
lar motallicitics (Pcntcr icci et al."^ 2002) and CO emission 
tBcrtpldi ct al. 2003) have been observed in the highest- 
redshift quasars at 2: ~ 6, dust is likely to be present 
in the ISM of even high-redshift galaxies. Because of 
their long scattering path-lengths, Lya p hotons are ex- 
tremely v ulnerable to dust attenuation jNeufeldl Il99d: 
ICharlot fc Fall ,1991.1 . and it was thought that this could 
account for low observed Lya equivalent widths com- 
pared to t hat expected from optical Balmcr emission lines 
iMeier fc Tcrlcvich 1981 : Hartmann ct al.,..198&L as well as 
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early failures to detect high-redshift galaxies in blank sky 
surveys. However, further work has shown that dust con- 
tent is not strongly correlated with Lya equivalent width 
(where dust content can be inferred from metallicity or sub- 
millimeter emission). For instance, some dust-rich galaxies 
have significantly higher Lya photon escap e fract ions than 
less dusty counterparts (| Kunth e t al.l|l99a . l2003ll . Indeed, 
iGiavalisco. Koratkar. fcCalzettTi^Off Tfound a lack of cor- 
relation between the equivalent width of Lya and the UV 
continuum slope l3, which measures continuum extinction. 
They interpreted this as evidence for decoupling of the ex- 
tinction of continuum and resonant line photons. 

Such decoupling could take pla ce if the interstel- 
lar medium is clumpy. iNeufeldl (Il99ll) and ICharlot fc Falll 
il993l) emphasized the importance of the geometry and 
multi-phase nature of t he ISM i n affe cting the observed 
Lya line. In particular, iNeufeld! (^9^) showed that in a 
clumpy, dusty ISM, the emergent Lya emission could have 
a higher equivalent width than the unprocessed spectrum 
of the underlying stellar population. For instance, if the 
dust survives primarily in cold neutral clouds, Lya pho- 
tons scatter off the clouds and spend most of their time in 
the intercloud medium, whereas continuum photons propa- 
gate unhindered into the clouds and suffer greater extinc- 
tion. Observationally, the ISM of our Galaxy is known to 
be clumpy down to small scales jMarschcr, Moore, & B aniai 
Il993l : IStutzki fc GuestenI Il990[l . with a power-law cloud 
mass spectrum b ased on GO (|g^n (^^rs. gc^ oville. fc SolomonI 
Il985l) and 21cm llDickev fc Garwood! 11989^ emission data. 
From IRAS lOO/xm, CO and 21cm data, there is evidence 
for a multi-scale fractal stru cture for both the diffuse HI 
clouds l |BazelLfc-^eserj|l98j) as well a s the molecular com- 
ponent ^El^^gree^^^^l^^on^ ^9^. The dumpiness of 
the ISM is well-established and it must be taken into account 
in radiative transfer calculations. 

Surprisingly, there has not been any detailed, quanti- 
tative, three-dimensional study of the effects of a dusty, 
cl umpy, ISM on L ya radiative transfer. The pioneering work 
of lNeufeldHl99ll) was a semi-analytic calculation for a plane- 
parallel slab: many issues, such as the detailed line pro- 
file and the effect of geomet ry, canno t be a ddressed with 
such an approach. Recently, iRichlin j i2003l) made a first 
attempt at a quantitative, three-dimensional calculation, 
but the slow convergence of the numerical technique em- 
ployed restricted the study to line center optical depths of 
r ^ 100, corresponding to neutral hydrogen column densi- 
ties oi N ^ 10^^ cm~^ for velocities v ~ 100 kms~^: orders of 
magnitude too low to be applicable to high-redshift galaxies. 
There have been many studies of the radiative transfer of U V 
continuum photon s in a clumpy, dusty ISM, using a variety 
of techniques Ce.g.lWitt fc GordonI (Il996l) : IVarsoi fc Dwekl 
il999l) ; iGordon et all i200ir) ). but none with extensions to 
resonance line photons. Conversely, while there have been 
Monte-Carlo radiat ive transfer studies of Lya photons in 
both static media jAhn. Lee, fc Led 200lL I2OO2,) and ex- 



panding supershells jAhn. Lee, fc Lejl2003l) . all have only 
considered a uniform medium. This paper therefore repre- 
sents a first attempt at numerically investigating Lya radia- 



tive transfer incorporating both the effects of dust and gas 
clumping. 

A key motivation is understanding recent puzzling ob- 
servations of anomalous equivalent widths in high-redshift 
galaxies. For instance, h i gh red shift 2 = 4.5, 5.7 sources ob- 
served by iRhoads et alJ i2003r) in the Large Area Lyman 
Alpha survey (LALA) show anomolously large Lya equiv- 
alent widths of EWJS 150 A (rest-frame), many far in ex- 
cess of any known nearby stellar population. An AGN ori- 
gin is unlikely, as the observed upper limit on the X-ray 
to Lya ratio is about 4-24 times lower than the ratio for 
known type II quasars {Malhotra & Rhoads 200 1)- The ra- 
diative transfer effects studied in this paper can produce an 
anomolously large Lya equivalent width from a standard 
stellar population. Such an effect could also be at work in 
the mysterious Lyman-alpha emitters observed at 2 ~ 3.1 
byiteidcl et al] feoOOt) . which have enormous Lya fluxes of 
~ 10~^^ergs~^cm~^ (a factor ~ 20 — 40 times large than 
typical line emitters at the same redshift), but no observed 
continuum. Finally, our calculation could be of particular in- 
terest in interpreting the large (^ 1 000) sample of Ly-break 
galaxy spectra (Shaplcv ct al1 l2003ll . as well as understand- 
ing the spectra of galactic starbursts with winds. 

The outline of this paper is as follows. In SJH we derive 
the basic multiphase Lya scaling relations. We then consider 
radiative transfer off opaque gas surfaces in ^ describing 
the Monte Carlo simulations and obtaining fitting formulae 
for the absorption probability, angular and frequency redis- 
tribution functions for both continuum and resonant scat- 
tering. With these surface scattering formulae in hand, we 
then develop a framework for multi-phase radiative transfer 
in [21 where we derive escape fractions and Lya line widths, 
discuss the role of the gas geometry, analyze several geome- 
tries of astrophysical relevance, and discuss the effects of 
dilute gas in between the opaque clumps. The surface scat- 
tering formulae substantially reduce the computational cost 
of simulating Lya transfer, making otherwise intractable cal- 
culations feasible. We also develop a simple analytic model 
with a single geometric parameter that shows good agree- 
ment with the full Monte-Carlo simulations. In ^ we dis- 
cuss some applications of our formalism. We show how pref- 
erential absorption of continuum photons can lead to strong 
enhancement of the Lya equivalent width. We also consider 
the typical Lya line profiles resulting from outflows/inflows 
of multiphase gas, and relate the profile characteristics to 
the outflow/inflow speed and the gas geometry. 



2 SCALING RELATIONS FOR Lya 
ABSORPTION 

In this section we build some physical intuition, by making 
simple order-of-magnitude estimates for the absorption of 
Lya photons in both homogeneous and multi-phase media, 
and summarizing some of the most important results from 
fj^and [2] We shall see that for conditions prevailing in most 
galaxies, Lya photons cannot escape unless the medium is 
multi-phase. 

Before beginning, it is useful to define some terms. Let 
fdop = (V^"^ /cjuo be the line Doppler width, where i/q is 
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Parameter 



Value 



HI line center resonant ctq 5.90 X 10 ^^T^ cm-^ 
scattering cross section 



Dust interaction cross a'^ 
section per hydrogen 
nucleus 



Dust absorption cross a" SdO"'* 

section per hydrogen ""-21 = "'"/lO^^"'^ cm^/H 

nucleus 



Absorption parameter f3 ct^/xhicto 



Damping parameter 



= 4.72 X 10~''T4"^/2 



Frequency in doppler x — vq) /^dop 

units 



Voigt function 



|x|>3 



Doppler speed 



ydop y^ii^ = I2.85T4I/2 km/s 



Absorption albedo 



Scattering asymmetry g (cosSgcat) 
parameter 



Table 1. Common Radiative Transfer Parameters Note: 
T4 = T/10^ K, Xhi is the hydrogen neutral fraction, vq = 2.48 X 
lO^^s"^ is the Lycf line center frequency, ^'dop = {V'^°'^/c)vq is 
the doppler frequency, i^l = 4.03 X W~^vq is the width of the 
quantum broadening Lorentzian profile, and Sscat is the angle 
between the incident and outgoing photon directions. 



the Lya line center frequency, and V'^°'^ = (2kBT/mp)^^^ 
is the characteristic atomic velocity dispersion times \/2- 
We evaluate the frequency shift from line center in Doppler 
units, X = (u — Vo^jvAop^ ■ The Lya scattering cross sec- 
tion is oix) — ao^{x), where is the Voigt function, 
which is characterized by a Gaussian Doppler core, and 
Lorentzian damping wings due to quantum broadening. For 
frequencies in the line wing, > 3, the Voigt function 
is dominated by the Lorentzian: "l>(a;) « a/(y'7Fx^), where 
a = !/L/2i^dop = 4.72 X 10~''^4~^/^ i^l = 4.03 x 10"*i^o is 
the width of the Lorenztian profile, and T4 = T/IO^K. For 
ease of reference, we have listed the most common radiative 
transfer parameters used in this paper in Tabled 



2.1 Homogeneous Slab 

We begin by reviewing the physics of radiative transfer of 
Lya photons through an o ptically thick sl ab, a problem that 
was first correctly solved bv lAdam!jJl972D . and subsequently 
verified and explored in much greater detail (Harrington 
19731 : iHummer fc KunasdllQSOl: iBonilha et alJli97fll: iFrisch I 



m;^. ■ IMummer fc KunaszlllijSd: Itjomllia et alJIlHYt* li^riscn I 
198(11: lNeufeldlll99nl: lAhn. Lee, fc Leell2002^ . We use these 



^ For gas at lO** K and a central frequency of 1216 A, 
the frequency and wavelength conversions are: 1 Doppler 
width=12.85 km/s=0.16A. 



classical results to test our Monte-Carlo code in Appendix 

El 



2.1.1 Homogenous Slab: Dust-Free 

Consider a Lya photon escaping from a dust-free slab of 
pure HI with line center optical depth tq. When the photon 
is the Doppler core, its mean free path is very short, and it 
barely diffuses spatially. It is always scattered by atoms with 
the same velocity along its direction of motion as the atom 
that emitted it. On rare occasions, it will encounter a fast 
moving atom in the tail of the Maxwellian velocity distri- 
bution, with large velocities perpendicular to the photon's 
direction. When this photon is re-emitted, it will be far from 
line-center, where the slab is optically thin. For a line-center 

optical depth of tq = 10^, a frequency shift oi x ~ 2.6 is suf- 

2 

ficient to render the slab optically thin, r ~ roe"^ ~ 1, 
and the photon can escape. So escape from the medium is 
dominated by rare scattering events. 

However, if the medium is sufficiently optically thick, 
Toa > 10^, the non- negligible optical depth due to the damp- 
ing wings still prevents escape. In this case, the photon will 
suffer repeated scatterings in the Lorentzian wings of typical 
atoms, and diffuse slowly in space and frequency, executing a 
random walk. Each scatter induces an r.m.s. Doppler shift of 
order x ^ 1, has a mean Doppler shift per scatter of — l/|x 
(with a bias to return to line center, due to the large proba- 
bility for photons to scatter in the core ; Ostcrbrock ( 1962)), 
and transverses an optical depth to${x)At ~ 1, or a mean 
free path which is At ~ l/"l>(a;) line-center optical-depths. 
Hence, a photon at frequency ja::j ^ 1 returns toward line 
center after jVe ~ x^ scatterings, having travelled an r.m.s. 
optical depth T^mB ~ v^A/TAt ~ \x\/^{x). If on its single 
longest excursion, the photon diffuses an r.m.s. distance of 
order the system size, Trms ~ |a;|/<E>(x) ~ tq, then the photon 
can escape. Since ${x) ~ a/x^, this implies a critical escape 
frequency: 

= (aro)^/^ « 30 T7l/'iV2V^ 

or almost ~ 400N2('^T4^^^ kms" 
where 7V21 = 7Vh//(10^^ cm~^). This displacement of pho- 
tons away from line center can be seen in our Monte-Carlo 
simulations in Fig. EH 



(1) 

away from line center. 



2.1.2 Homogeneous Slab: Dusty 

Now let the gas contain dust, with an total (scattering -f ab- 
sorption) interaction cross section per hydrogen atom of cr'*, 
and an absorption probability per dust interaction td- The 
average absorption probability per interaction with either 
dust or hydrogen is: 



^absorb 
""total 



Xiii${x)ao +a'^ 

-3 T4 (7^.21 



(2) 



<i>{x) 



1.59 X 10" 



Xhi 



where Xhi is the hydrogen neutral (HI) fraction (which must 
be introduced because a'' is the cross-section per hydrogen 
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nuclei). In the third step we used the fact that, except very 
far from hne center, HI scattering dominates, Xiii^{x)ao 2> 
a'', and defined the absorption parameter 



/3 = €da/xHicro = 1.69 x 10 



-8 



^HI 



(3) 



where o"" 
Way, al-)-, 



td(J ■ For the diffuse HI phase of the Milky 



ff VlO" 



IHja 1, td ~ 0.5, and Xm « 1 



JPraine fc Lee Ill984 Iwhittet Ibooi iDraine 112001 . and so 

f} ~ 10~*. The fourth step in eqn. l(2j uses the wing photon 
approximation ^(x) ~ a/^/TTx^. 

Under what conditions can the Lya photon escape from 
such a dusty medium? While this has be en the subject of de- 
tailed a nalytic and numerical work (e.g ., lHummer fc KunasJ 
Jl980l) : lFriscii1 il980l) : iNeufeldl Jl990l) '). we can understand 
the basic scaling laws quite easily. The probability that a 
photon will be absorbed at a given frequency x is simply 
the number of scatterings at that frequency times the prob- 
ability of absorption per scattering: 



Pibf ~ J^{x)e{x) 
Thus, Plifix) ~ 1 for 



X p — . 

a 



\X\ > Xs,bs 



1/4 



12.9 



Xhi 



1 1/4 



T4 (Tt- 



(4) 



(5) 



This implies that a photon will be absorbed before es- 
cape if it has to diffuse far into the line wings in order to 
escape from the slab, or if 

Xc ^ 3^abs, 

with Xc and a;abs given by eqns. Q and respectively'^. 
Hence, if the line center optical depth exceeds a critical 
value. 



To > Tc 



a/33 



1/4 



4.6 X 10« T-'^^[x^i/aUi?^\ (6) 



photons cannot escape from the medium. This simple cri- 
terion is borne out by mor e deta iled calculations (e.g., see 
Fig. 5 of lAhn. Lee, fc Lee I ll200Cfl . and references therein). 
In terms of the HI column density A'21, Lya photons cannot 
escape from a homogeneous dusty slab once 



N21 >o.08r4''"[xHi/(7^2i] 



13/4 



(7) 



Since typical HI column densities in the Milky Way and 
other galaxies is ~ 1, Lya photons could not es- 

cape if most of the HI is in a homogeneous dusty slab. 
In the next section, we see that if the gas is instead 
inhomogeneous/multi-phase, Lya photons can escape much 
more easily. 



^ Note that at Xabsi the absorption probability per interaction 
is still small, e <C 1, and scatterings still strongly predominate. 
The scattering and absorption cross-sections are comparable only 
at a much larger frequency, x(e ~ 0.5) ~ (a/P)^^'^ ^ S> 
aJabsorb, by which time all photons have been absorbed. 



2.2 Multiphase Gas 

We shall now estimate the absorption criteria for a multi- 
phase dusty HI distribution. It is worth first noting that a 
medium is always more transparent when it is clumpy, for 
fairly generic and model independent reasons. The effective 
optical depth in an inhomogeneous medium is rdumpy = 
—In ({exp [— t])), where the average is over all lines of sight. 
However, for a uniform medium, r = constant along all lines 
of sight, so that Tunitbrm = (t) . From the standard triangle 
inequality. 



{exp[-r]) ^ exp[-{r)] 



(8) 



and applying the negative logarithm to both sides, we see 
that Tciumpy ^ Tuniform . Thus, for iustauce, flux transmission 
in quasar absorption spectra is increased for an inhomoge- 
neous intergalactic medium (IGM), where transmissi on is 
dominated by underdense voids (e.g.. iFan et al I (|2002f )'l. 

This effect is strongly exacerbated if most of the ab- 
sorbing material lies in dense clumps which are optically 
thick to scattering. In this case, most of the photons scat- 
ter off the cloud surfaces without penetrating the clouds, 
which effectively shields the absorbing material. This situa- 
tion naturally arises in a multi-phase ISM, when most of the 
dust lies in dense molecular/atomic clouds. For now, let us 
assume that the inter-cloud medium is highly ionized and 
relatively dust-free, so that all of the dust and HI lies in 
dense clouds. 

In H3.3I we show that we can calculate analytically the 
escape probability of Lya photons in a multi-phase medium 
quite accurately, given just two parameters: Ao and e^. We 
define Ao as the mean number of cloud surfaces a photon 
would encounter before escape in the absence of absorption. 
It only depends on the geometry of the multi-phase medium 
(the trajectory of photons is independent of frequency, pro- 
vided clouds are very optically thick). For most of the cases 
we will consider, Ao ~ 1 — 30, with typical values Ao ~ 5. 
The cloud albedo Ec is the probability of absorption upon 
hitting a cloud surface. In H3.3I we show that: 



(9) 



where e is the absorption probability per interaction given 
by eqn. ((^J. Eqn. Q is easily understood in the case where 
e is constant (e.g., for coherent scattering). The effective 
absorption optical depth of a medium with scattering is 
r, ^ \/T^r(Ta + T^ (e.g., Rybicki fc Lightman, 1979, p. 38), 
where and are the absorption and scattering optical 
depths, respectively. Hence, the albedo is tc = t*/ (rs+Ta) ~ 
\/ Tajija -f Ts) = y^. lu H3.3I wc show this scaling still holds 
for Lya photons, despite the fact that eix) changes as the 
photon random walks in frequency whilst scattering within 
the cloud. We find that the typical frequency shift after scat- 
tering off an optically thick surface is Ax ~ 1.5, with most of 
the redistribution in a symmetric profile about the incident 
frequency Xi. With e(j:) given by eqn. (|5J, the symmetric fre- 
quency distribution about Xi implies ^ t{x)^ ~ \/e{(x)) ~ 

^ e{xi), since (x) = Xi. Therefore the coherent scattering 
absorption law describes the Lya absorption when e is eval- 
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uated at the incident frequency, ~ 2^/ e{xi). We verify 
this explicitly in t|3.3.2l 

Under the above approximations, the probability that 
a Lya photon at frequency x will be absorbed is; 

1/2 

|a;| (10) 



pmulti-phasc _ j^^^^ _ 2AAo 



which should be compared against eqn. I^J for ^ homoge- 
nous slab. Notice the much weaker scaling with frequency: 



rimulti — phase 
abs 



oc X, instead of PfJ,? 
P'^''^'^ ~ 1 for frequencies 



From eqn. HUH . 



1/2 



(11) 



16.7 



Xhi 



1/2 



T4aUi. ' ^^^^ 

Comparing against eqn. JSJ, the cut-off absorption frequen- 
cies for the slab and multi-phase case are actually compara- 
ble, modulo the value of Ao- However, there is an important 
difference: Lya photons have to diffuse far into the wings in 
order to diffuse spatially out of an optically thick slab. Since 
escape requires x^ > a;absorb for a very optically thick slab, 
the photons will inevitably be absorbed. By contrast, there 
is generally much less diffusion into the line wings when 
scattering off surfaces in a multi-phase medium. Photons 
typically only penetrate small optical depths, r ~ 1 — 10 in 
the cloud surfaces before escaping, and the number of scat- 
terings is much less. Thus, the majority of photons need not 
necessarily stray far from line center. 

Clearly, the crucial parameter which determines if Lya 
photons can escape from a multi-phase medium is Xe, the 
characteristic escape frequency; this must be small, Xc < 
i^absorb, for photous to cscape. Lya photons in a multi- 
phase medium acquire Doppler frequency shifts in two ways: 
through the thermal motions of HI atoms, as before, and 
also through the bulk motions of the clouds/scattering sur- 
faces. For this reason, it is useful to rewrite a^abs in units of 
velocity: 



Vo""^ = 2.1 



[xai/al 



,1/2 



(13) 



where V2 = 1//100kms~\ In SjO] we show that atomic 
motions cause a net r.m.s. frequency shift of ~ l.SAo 1^''°'' 
after Ao surface scatterings, and consequently result in an 
escape frequency of: 



O-bV^^J^o ~ 0.3 [Ti]^^^ [Afo/5] 



(14) 



By contrast, we find that cloud motions (either random mo- 
tions or bulk inflows/outflows) with characteristic velocity 
causes a net r.m.s. frequency shift of: 

,(^^e,cloud _ _ 2.2 V2 [JVo/5\^^^ . (15) 

Note that y<=>="°-"i<^ oc A^o, while V<=><^i°"ds ^ ^yjj^^ ^j^j^j^ 
discuss in H4.3I When both the frequency redistribution and 
surface motion are combined, we find that the typical escape 
velocity is simply given by the sum (rather than the sum in 
quadrature). 



(16) 



For absorption to be important, we require V > 
ij^absorb^ which coustraius At) to be 

A-o>5 [0.38{T4 al2i/xui)'^^ + {[V2'fcTl2i/xni)'^Y' '(^^^ 

If Ao satisfies this inequality then the Lya photons will be 
significantly absorbed. This approximate constraint has the 
correct limits when either V"'"'^"'^''' = or V"'"^"""^ = 0, but 
is off by a factor of ~ 2 when V"'"'^""'''' ~ 't/<=><^i°"d -phe geo- 
metric parameter A/o thus plays a key role in determining if 
Lya photons can escape in a multi-phase medium, and plays 
an analogous role to the column density A^hi in a homoge- 
neous slab. In H4.2I we provide formulas for Ao for various 
different basic types of multi-phase geometries. 



3 SURFACE SCATTERING & ABSORPTION 

Multi-phase radiative transfer typically involves photon 
propagation through an optically thin inter-cloud medium, 
and repeated scattering off optically thick clouds. In a full- 
blown Monte-Carlo simulation, the latter consumes by far 
the lion's share of computational time. This is extremely in- 
efficient: the same scattering/absorption problem off cloud 
surfaces is being solved over and over again for each photon. 
A better approach is to consider each cloud as a scatter- 
ing/absorbing particle with its own radiative trans fer prop- 
erties (for other applications of this viewp oint, seelNeufeldl 
( 19911 : .Hobson fc Padman (,.1993.' l: .Varsoi fc Dwek Hl99gl) '). 
Wo characterize these cloud scattering properties in this sec- 
tion. 

For Lya photons, clouds are extremely optically thick 
and have essentially the same radiative transfer properties 
as a semi-infinite slab. This eliminates detailed dependence 
on the geometry of the cloud: all that matters is its dust 
content and the initial photon frequency. Surprisingly, the 
radiative transfer properties of a dusty semi-infinite slab to 
Lya photon scattering has not been characterized in de- 
tail. We do so in this section. We derive formulas for the 
net absorption probability (the "cloud albedo" ) Ec , the exit- 
ing photon angular distribution D(9), and the exiting pho- 
ton frequency redistribution R{xi,x), as a function of the 
initial photon frequency Xi and the gas composition. With 
these surface transfer formulae, radiative transfer through 
regions containing opaque gas clouds can be quickly esti- 
mated and/or simulated without performing any scattering 
calculations within the individual gas clouds. This allows for 
both vast speed-ups of Monte-Carlo simulations (outlined in 
>I4.61 and a tractable analytic multiphase radiative transfer 
analysis (S2J- H ^ photon typically scatters Af times before 
exiting a cloud (where Af ~ lO^"'^ for incident frequencies 
Xi ~ 5 and a"2i ~ 1 — see H3.2.4II . then this allows a speed- 
up of order ~ A/", making tractable multi-phase calculations 
which would otherwise be prohibitively expensive. 

Our approach is to find fitting formulas to Monte-Carlo 
simulations of an ensemble of incident photons. Whenever 
possible, we base the fits on known analytic formulas for 
simple cases, extending the analytic formulas to encompass 
the more general cases that we simulate. We begin by de- 
scribing the Monte-Carlo algorithm we use. Surface radiative 
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transfer of continuum photons, where scattering by dust is 
effectively coherent, is discussed next. We then consider the 
more complex case of Lya surface scattering, where resonant 
frequency redistribution effects must be dealt with. Lastly, 
we discuss two kinematic aspects of surface scattering: the 
average scattering angle, and the frequency shift due to a 
bulk surface velocity. 



3.1 Monte-Carlo Code Description 

The Monte C arlo algorithm w e use is sim i lar t o the 
code used by 'Ahn, Lcc, & Lee| J200ll l200l l2003h and 
[Zheng fc M iralda-Escudc ( 2002J. These papers provide a 
fuller description of the algorithm than that which we give 
here. An ensemble of photons is run through a medium with 
neutral HI and dust, and statistics are gathered. Each pho- 
ton is tracked until it either escapes the medium or is ab- 
sorbed, at which point the photon's flight is terminated. The 
optical depth r between each interaction is drawn from the 
distribution exp(— r), i.e. r — — Inu where u £ [0,1] is a 
random variable drawn from a uniform distribution (here- 
after "univariate"), and the photon's position is updated. 
The photon then interacts with either the HI or the dust, 
resulting in either HI resonant scattering, dust scattering, 
or dust absorption, all of which we describe next. 

We model the dust as particles which can either absorb 
or coherently scatter photons. Although dust scattering is 
not necessarily coherent, in practice ignoring frequency re- 
distribution due to dust is an excellent approximation. The 
trajectory of a continuum photon is unaffected by small de- 
viations from coherent scattering, since the dust albedo ed 
only varies weakly with frequency. However, the Lya ab- 
sorption probability is a strong function of frequency (see 
eqn. (|5J, so for resonant scattering the effects of frequency 
redistribution must be taken carefully into account. 

A continuum photon only interacts with dust, with an 
absorption probability ed per interaction. We determine if 
the photon is absorbed during a given interaction by draw- 
ing a random variable it e [0, 1] from a uniform distribution; 
if It ^ ed, the photon is deemed to be absorbed. If the con- 
tinuum photon is not absorbed, then its direction is changed 
by the scattering angle Sscat off the incident direction, with 
a random azimuthal angle. We use the H enyey-Greenstein 
scattering angle distribution^ (e.g., dWittl l 19771') 



;(e; gd) 



9d 



An 



[l + gd - 2gd cos^ 



-3/2 



(18) 



which is parameterized by the dust scattering asymmetry 
parameter Qd € [—1,1], and where we use the normaliza- 
tion 1 = Jq dd P{0). The scattering asymmetry parameter 
is defined as g = (cos^scat). To appr oximate dust absorp- 
tion and scattering in the Milky Wav iDraine fc Lee 11198^: 



^ iDraine 1 12003^ shows that the Henyey-Greenstain distribution 
is inaccurate for wavelengths A ^ 4700 A. However, since the sur- 
face scattering problem we consider is essentially planar, the de- 
tails of the dust scattering distribution should significantly affect 
the results. 



IWitt fc GordonllOQdIwhittet l2003l:lDraine l2003l) at wave- 
lengths near 1216 A, we use = 0.5 and gd = 0.5, unless 
otherwise noted. 

If the photon is a Lya photon, then the interaction can 
either be with dust or with neutral HI. The probability of a 
dust interaction is a''' /{a'''-\-Xiii^{x)ao); if a random univari- 
ate u is less than this, then the interaction is identical to the 
"continuum" dust interaction described above. Otherwise, 
the Lya photon scatters resonantly off neutral hydrogen. In 
all our simulations we take the hydrogen in the cold phase to 
be completely neutral, and so adopt a^ni = 1 unless otherwise 
noted. Although the velocity distribution of the hydrogen is 
Maxwellian, the velocity distribution of atoms that scatter 
photons depends upon the frequency of the photon. Let h 
be the direction of the photon before scattering. In the two 
directions perpendicular to h, the scattering atom's velocity 
distribution is Gaussian, 



/(wi 



(19) 



where w± = v±/V'^°'^ and v± is a velocity component in 
one of the two transverse directions to h. In the direction 
parallel to n the velocity distribution of scattering atoms is 
a Gaussian weighted by a Lorentzian: the Gaussian is due to 
the thermal motion of the gas, while the Lorentzian is due 
to the increased probability for scattering in the wings from 
quantum mechanical broadening. The velocity of the atom 
Vz along the direction of the incident photon is determined 
by drawing a random variable from the distribution: 



a exp(— tt)f) 



where = v^/V'^°'^ fsee lZheng fc Miralda-Escudel JioO 



(20) 



for a rapid algorithm for generating random numbers with 
this distribution) and Xi is the photon's incident frequency. 
In the rest frame of the atom, the frequency of the outgoing 
photon is the same as the incident frequenc y (strictly sp eak- 
ing it differs slightly due to the recoil effect ^Fieldll959^ , but 
for our purposes this is negligible). The new direction h' is 
given by a dipole distribution, with the symmetry axis de- 
fined by the incident direction n: 

p{e) = ^ii + cos^e) , (21) 

where 9 is the polar angle off the direction h. Although 
resonant scattering can result in either isotropic or dipole 
scattering angle distribution s, depending upon the interme- 
diate excited quantum state JStenfidll980l:lAhn. Lee, fc Led 
l2002h . the difference is immaterial for calculating spectra 
and escape fractions; a more careful treatment would be re- 
quired, for instance, to accurately simulate Lya polarization. 
Given the new photon direction h' and the scattering atom 
Doppler velocity w, the new photon frequency x' is given by 



h ■ w + h ■ w 



3.1.1 Avoiding Core Scatters 



(22) 



Lya photons spend most of their scatters in the line core, 
where spatial diffusion is typically negligible. Essentially, 
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each time a photon enters the hne core it scatters in place 
until it is scattered by a high speed atom which moves the 
frequency out of the core. The frequency at the core-wing 
boundary, Xc ~ 3, is defined by: 



1 



(23) 



Hence, it typically takes exp(a;c^) ~ exp(9) ~ 10* scatters 
to scatter out of the core. Note that Xc depends only log- 
arithmically on the gas temperature, through the damping 
parameter a. For a photon that starts out at frequency Xi 
in the line wing, the photon typically returns to the core 
fairly quickly, after ~ Xi'^ = 25{xi/5)^ scatters (see discus- 
sion in tl2.H . Consequently, most of the simulation time is 
spent calculating core scatters. By circumventing the core 
scatters, the simulation can be greatly sped up. We have 
adopted a scheme to d o so that is similar to that used by 
lAhn. Lee, fc Le3 i2002h , with the addition that we also con- 
sider absorption. 

One might think that absorption whilst scattering in 
the line core is negligible, due to the small physical path 
lengths traversed whilst scattering in the core. We confirm 
this quantitatively below. 

For a photon with an initial core frequency Xi where 
|a;i| < Xc, let fie be the average number of scatters for the 
photon to leave the core, Xc = be the average frequency 
(absolute magnitude) while in the core, and x^ = {\x-w\} 
be the average frequency (absolute magnitude) of the first 
scatter that leaves the line core, x„ > Xc- We ran simulations 
for gas at 10* K and adopted Xc = 3. We find that for any 
initial frequency in the core, fie ~ 2.9 X 10*, Xc ~ 0.57, and 
Xw « 3.3 with an equal probability of leaving the core at 
X = 3.3 and x = —3.3. Thus, the probability of absorption 
during core scatters can be approximated by (see eqn. @): 



4.0 X 10*/3 = 6.8 X 10"* [xh 



(24) 



0"-21 



Since the typical photon scatters A/" ~ 10 times before es- 
caping the surface f H3.3.4l l and it takes Af""'" ~ 9{xi/3f 
scatters to reach the core, a typical photon injected in the 
line wing will visit the core perhaps once. The probability 
of a Lyof photon being absorbed in the core during the sur- 
face scattering is, therefore, ~ 10~*(t" 21 (^d/CS), which is 
negligible. 

We therefore devised the following acceleration 
scheme*: 1) If a photon starts off in the line core, we do the 
exact core scattering, and only employ the approximation 
scheme on subsequent visits to the core. This is computa- 
tionally cheap, since an incident core photon does not pen- 
etrate deep into the surface, and typically leaves after a few 
scatters. 2) If a photon enters the line core from the wing, 
the probability of absorption is ecorc. 3) If the photon is not 
absorbed, then it is given a wing frequency x = ±Xjn, with 
an equal probability for plus and minus. The spatial position 

* For completeness, this scheme still takes core absorption into 
account. This may be useful in other contexts when the core is re- 
visited many times and core absorption could be non-negligible — 
e.g., for Lya photons escaping from an optically thick slab. 



of the photon is exactly the same as where it entered the line 
core, and the new angular direction is randomly drawn from 
an isotropic distribution. In i|A2l we compare this acceler- 
ated scheme to exact simulations of surface scattering. In 
practice, it gives accurate results, and gives a vast speed-up 
of the simulations, typically of order ~ 10^. 



3.2 Surface Scattering of Continuum Photons 

In this section, we study the properties of coherent surface 
scattering of continuum photons. It is very useful to under- 
stand the properties of coherent scattering surface transfer 
in order to have a baseline for comparison with Lya surface 
transfer. As such, in this section we do not use the Henyey- 
Greenstein scattering angle distribution, eqn. II8II . but in- 
stead use the same distribution that we use for Lya scatter- 
ing, which is the dipole distribution, eqn. 1211 ^. We begin 
by calculating how thick a slab of gas must be before the 
surface scattering approximations apply. We then describe 
fits for the cloud albedo ec, the exiting photon scattering 
angle distribution D(6), and the typical number of scatters 
N. 



3.2.1 The Surface Approximation 

For a slab of material with a finite optical thickness r, the 
radiative transfer of photons incident on a surface will be 
approximately the same as for a semi-infinite slab if the frac- 
tion of photons that are transmitted through the slab, 
is small. In this limit, the surface is not translucent but acts 
as an absorbing mirror, and all photons are either reflected 
or absorbed. We define the penetration column density A''* 
such that when N > N^^ the transmitted fraction is less 
than 10%, < O.l. From a series of Monte-Carlo simula- 
tions, we find that a decent fitting formula for /* is 

f = [(1 + r) cosh (^0.55r'/*e'/2)] , (25) 

as shown in Figure Q The transmission will be negligible 
(/^ ^ 0.1) when either re^^'^ 3 or r ^ 9. The corresponding 
penetration column density is 

< = "ii'if-T^. -^"l • (26) 



3.2.2 Surface Absorption 

To derive a formula for the cloud albedo ec, we ran simu- 
lations for an isotropic surface source and averaged the ab- 
sorption over this ensemble. For one dimensional radiative 
transfer, an exact formula for ec can be derived for photons 
incident on a semi-infinite line of material, 



(27) 



^ Note that when we actually perform Monte-Carlo simulations 
of Lya photons, we do use the Henyey-Greenstein distribution 
when the Lya photon scatters off dust. 
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Figure 1. Transmission of Continuum Photons The fraction 
of incident photons that arc transmitted through a finite slab, 
is shown as a function of the slab's total optical depth r, 
for various values of e. The circles are simulation results; from 
lightest to darkest, e = 0, 0.01, 0.1, 0.5. The thin lines are the 
fitting formula eqn. 1251 . The thick line is the limiting case = 
exp(— r), which corresponds to e = 1. 




Figure 2. Continuum Photon Surface Absorption Surface 
absorption probability Ec for continuum photons, as a function of 
dust albedo e. The diamonds are simulations for coherent scat- 
tering with a dipole scattering angle distribution. The line is the 
fit eqn. (Til . 



where each scatter is front-back symmetric {g = 0) . As 
shown by Figure|5| eqn. 1271 provides a very good fit for the 
three dimensional, semi-infinite plane case. When e ^ 1, we 
find ~ 2^/e . which is similar to the power law found by 
lNeufeldHl99ll) . 

3.2.3 Escape Angles 

To find a fit for the distribution of exiting angles for photons 
that escape, we ran simulations for various incident angles 8i 
relative to the surface normal. We find that for nearly all 6i, 
and for both isotropic and dipole single particle scattering, 
the distribution of exiting angles is well fit by 

Dss{0) =sm2e, (28) 

as shown in Figure |3 This distribution can be understood 
as the combination of two effects. First, if photons scatter 
multiple times before escape, the photons effectively lose all 
"memory" of the incident angle. This effect leads to a ran- 
dom exiting angle distribution, Dss{d) oc sin^. Second, pho- 
tons that exit with angle 9 will be attenuated if the optical 
depth traversed during the exiting leg, r, exceeds unity. Let 
be the perpendicular optical depth at the point of last 
scatter for a photon that would exit with an angle 9 in the 
absence of absorption. The condition r 1 implies a max- 
imum perpendicular depth of such a photon is r^f^ax = cos 9 
(e.g., only shallow surface layers contribute photons escaping 
nearly parallel to the surface). Near the surface, the mean 

^ For a plane-parallel slab, the Eddington approximation with the 
two-stream boundary condit ion gives /e oc \/e/{l + (see, e.g., 
jRvbicki fc Lightman ][l973), pg. 320), and from our simulations 
we find that the prefactor is 2 for 1-D scattering. 



intensity is approximately constant for r <, 1. This implies 
that the number of photons available to escape at an angle 
9 scales as cx Tj^axi which implies Dss{9) oc cos^. Includ- 
ing both the effect of randomization and attenuation gives 
a distribution Dss{9) oc sin cos 6', which, when normalized 
over 9 £ [0, 7r/2], gives eqn. I|28^ . In the case of dipole scat- 
tering there are deviations from this fit for grazing incident 
angles, but overall this fit is generic for surface scattering 
when the single particle scattering distribution is front-back 
symmetric. We also find that there is very little dependence 
on the absorption albedo e. When the distribution Dss{9) 
holds, the distribution of azimuthal angles is uniform over 
the interval [0, 2n). Finally, we note that if the slab is viewed 
at an angle ^obs (from the outward surface normal), then the 
observed intensity is oc Dss{9ohs) cos Sobs. The extra cos Sobs 
factor is due to the dependence of the projected surface area 
on the viewing angle. 



3.2.4 Number of Scatters for Escape 

In H4.1l we present a general derivation of the average number 
of scatters. A/", as a function of the escape fraction fa and 
the absorption albedo e, given by eqn. I|55^ . Application of 
this formula to surface scattering, where the escape fraction 
is /o = 1 - ec(e), gives 

AA = , (29) 

which is shown by the solid line in Figure 2] As the ab- 
sorption albedo goes to zero, the average number of scatters 
diverges, although the median number of scatters does not 
seem to increase beyond ~ 5. Clearly, the average is domi- 
nated by the rare photons that wander deep into the surface. 
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Figure 3. Exiting Angle Distribution The distribution of 
exiting angles relative to the surface normal is shown for sev- 
eral incident angles 9i and absorption albedos e. The shaded re- 
gion denotes the distribution Dss{S) = sin 29. The thick lines 
are dipole scattering: from darkest to lightest, the lines depict 
(e,6»i) = (10-*, 0°), (10-*, 45°), (0.1,45°). The thin line is for 
isotropic scattering, {e,di) = (10—^,45°). 




10 10^ 

e , e(x.) 



Figure 4. Number of Scatters: Continuum & Lya Pho- 
tons The number of scatters for escape for coherently scattered 
photons (circles) and Lya photons (squares and triangles) as a 
function of the absorption albedo. For Lya photons the absorp- 
tion albedo is evaluated at the incident frequency; two initial fre- 
quencies are shown, Xi = 10 (squares) and Xi = 20 (triangles). 
For each case, both the average (filled symbols) and median (open 
symbols) number of scatters are shown. The analytic estimate for 
the number of scatters for coherent scattering (solid line), given 
by eqn. 1291 . matches the data very well. 



tribution D{9), the typical number of scatters, and the Lya 
frequency redistribution R{xi,x) as a function of the inci- 
dent frequency Xi. 



3.3.1 The Lya Surface Approximation 

The criterion for the surface scattering approximation to 
apply for incident Lya photons is similar to that defined 
for coherent scattering. Consider Lya photons with fre- 
quency Xi incident on a finite slab with optical thickness 
Ti = To^{xi), evaluated at the incident frequency. When Xj 
is in t he Lorentzian wing and if the slab is pure HL lNeufeldl 
ilQQOl) analytically derived that the fraction transmitted is: 



/hi — — 
in 



(30) 



when Ti 3> 1. Our simulations confirm this result when there 
is very little dust, but find that can be substantially less 
than /hj when a small amount of absorbing dust is present. 
To derive a fitting formula for we ran simulations for var- 
ious incident frequencies and dust absorption cross-sections. 
To a large degree, depends only upon the incident single 
scattering albedo and the incident slab optical depth r^; 
almost all the frequency dependence is captured by these 
two parameters. This is extremely convenient: the transmit- 
ted fraction (and as we shall subsequently see, the refiected 
fraction) depends in a fairly simple way on the properties of 
the slab and the incident frequency. One might worry that 
due to frequency redistribution (which can be substantial; 
see >I3.3.5L the frequency dependence becomes extremely 
complicated, but that appears not to be the case. In partic- 
ular, the typical incident photon never "loses memory" of 
its initial frequency. Most photons that escape do so after 
a handful of scatters, ~ 10. When Xi <^ 4, the majority of 
photons do not wander into the core before escaping, and so 
most photons retain some memory of the incident frequency. 
As shown in Figure |S] a reasonable fit for photons initially 
in the line wing is: 



3ri 



cosh ( [Tiy/nf'^'^'^ 



(31) 



When ei —> 0, the pure HI formula /hi, eqn. (1251 . is recov- 
ered when Ti ^ 1. As can be seen in the figure, for pho- 
tons initially in the Doppler core the transmitted fraction 
is slightly larger than this. The transmission will be neg- 
ligible if either riy^>,3 or >, 12. This corresponds to a 
penetration column density 

NF,l = min (o.2 VilaUiV'^^ 0.05 [Fa]') (32) 

for any incident frequency \xi\ ^ 3, with substantially 
smaller for photons in the Doppler core. 



3.3 Surface Lya Transfer 

As in the case of coherent scattering considered above, we 
begin by calculating how thick a finite slab of gas must be 
before the surface scattering approximations apply. We then 
describe fits for the net absorption ec, scattering angle dis- 



3.3.2 Surface Lya Absorption 

To derive a fit for the surface albedo Cc, we ran again sim- 
ulations for a wide range of incident frequencies and dust 
cross-sections. As in the coherent scattering case, we aver- 
age Ec over an isotropic incident direction. We find that Cc 
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Figure 5. Transmission of Lyo Photons The fraction of es- 
caping photons that are transmitted through the slab as a 
function of = e(xi), for gas at T = 10^ K with dust parameters 
Qd = 0.5 and = 0.5. Each shade is a different incident slab opti- 
cal depth Ti: from lightest to darkest = 1, 2, 4, 8, 16, 32. Sim- 
ulation data for three different incident frequencies shown: 
the diamonds, squares, and circles are Xi = 1, 5, 10, respectively. 
For each Xi and Ti, simulations were run for three different dust 
absorption cross-sections: reading from left to right for each Xi 
and Ti, the dust values are ~ ^-^^ 10- T^^^ lines are gen- 

erated from the fitting formula eqn. lijll . 



mainly depends upon the single scattering albedo at the in- 
cident frequency, ti. As shown by the solid line in Figure |S1 
Ec is well fit by 



(33) 



which has the same form as for coherently scattered photons, 
eqn. 12711 . A slightly better fit is shown by the dashed line 
in Figure |S] 



l + 2e,i/2 



(34) 



Since eqn. H34|l gives a slightly better fit only when tc is 
negligibly small, in practice we always use eqn. I|33|l : this 
will simplify the subsequent multiphase analysis somewhat, 
for only a small loss in accuracy. Either of these formulas 
for ec applies even when Xi is in the Doppler core (although 
eqn. 13411 gives the better fit in this case). 

As long as the Lya photon is in the line wing, the ab- 
sorption probability is independent of the gas temperature. 
To show this, in Figure |7] we calculate ec using eqn. 13311 . 
as a function of the incident velocity shift AV in physical 
units, rather than Doppler units x. The figure shows the cal- 
culation for gas at temperatures T — 100 K and T — 10* K 
and for varying dust content cr" 21- For photons in the line 
wing, the temperature dependence of e drops out: 



■<■■■' o 




XHi$(a;)cro ax-taOQ axuicro 



(35) 
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Figure 6. Lya Surface Absorption (1) The surface albedo 
Cc is plotted against the single scattering albedo evaluated at 
the incident frequency, , for gas at T = 10* K with dust pa- 
rameters gd = 0.5 and = 0.5. The inset shows an enlarge- 
ment of the e = 0.01 — 1 region. Simulations were run for four 
different dust absorption cross-sections and seven different inci- 
dent frequencies. Each shade corresponds to a different cr": from 
lightest to darkest cr^ji = 0.01, 0.1, 1, 10. For each value of cr", 
seven different incident frequencies Xi were run: from right to left, 
Xi = 20, 10, 5, 4, 2, 1, 0. Note that for Xi = 0, the value of for 
= 0.01 and 0-^21 = 0.1 is less than 10~^, and hence lies 
off the plot. The absorption probability Cc depends only on e^; at 
fixed ei, there is no independent variation with dust content cr'^21 
or incident frequency Xi. The solid line is the fit eqn. 1331 . while 
the dashed line is the fit eqn. 1341 . 



Focusing just on the temperature, a oc l/y/T, ao oc 1/\/T, 
and x = V/V'^°'^ cx l/VT. Therefore the combination 
x'^/aao has no dependence on the gas temperature. Physi- 
cally, scattering in the Lorenztian wing is dominated by the 
quantum broadening of the cross-section, which does not 
depend upon the thermal motion of the scattering atoms. 
Since we have shown that ec mainly depends upon ei, it 
follows that ec is also temperature independent. Thus, the 
surface absorption probability does not depend upon the gas 
temperature 

3.3.3 Escape Angles 

As argued in H3.2.3l the scattering distribution Dss{d), eqn. 

is fairly generic when the single particle scattering is 
front-back symmetric {g = 0). Since Lya scattering — which 
is either isotropic or dipole — is always front-back symmetric, 
the Lya surface scattering angle distribution is also given by 
Dss{d), independent of the incident angle 9i. 

3.3.4 Number of Scatters for Escape 

The average number of scatters for Lya photons to escape 
is more complicated than for continuum photons mainly be- 
cause Lya can be trapped in the Doppler core. As discussed 
in il3.1.1l Lya photons in the Doppler core must scatter 
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Figure 7. Lya Surface Absorption (2) Simulation results 
for the net probability of surface absorption plotted as a func- 
tion of the incident Lya frequency shift off line center, in veloc- 
ity units. The solid lines are for a gas temperature T = 10'' K, 
while the dashed lines are for T = 10^ K. We also show the ef- 
fect of varying the dust content: from the thickest line to the 
thinnest, ct^ji = 100' 10' 1' 0-1' 0-01, 0.001. The cloud albedo 
Ec is independent of gas temperature if the photon is in the line 
wing, though the core-wing boundary is obviously temperature- 
dependent: Vc ~ 'iV'^°^ oc T^/^ (the vertical dashed lines show 
the core-wing boundary at T = 100 K, 10** K). The lines flatten 
out once the incident frequency shift lies within the Doppler core, 
because the Lya scattering cross section approaches its line center 
value (TO • This causes the single scattering albedo to asymptote as 
y — > 0. Since <x ^Tt approximately holds even for incident pho- 
tons in the Doppler core, the surface albedo will also asymptote 
as y ^ 0. 



~ 3 X 10* times before a rare scattering event brings the 
frequency into the line wing. In Figure 01 we show exact 
Monte-Carlo simulations of Lya photons, where the core 
scatters are directly calculated, i.e. the Monte-Carlo accel- 
eration scheme described in is not used. The huge 
increase in scatterings over the continuum scattering result 
is obviously due to scattering in the Doppler core. The aver- 
age number of scatterings in the line wing is (analogous to 
the continuum scattering result) 



AA" ~ l/^/e{x,) ~ 25 [a;Hi/<7^2i]'^ [xr/^ 



(36) 



while the average number of scatterings required to reach 
the Doppler core is ~ 25[2;i/5]^ (see the discussion in H2.2II . 
Thus, although the typical photon does not reach the core, 
the probability of reaching the core is large enough that the 
the average N is dominated by core scattering. 



3.3.5 Surface Lya Frequency Redistribution 

In this section, we find a formula for the refiected frequency 
distribution R{xi,x) as a function of the incident frequency 
Xi, for the pure dust-free HI case. When dust is added, 
we find that the distribution i?dust(3::) adheres closely to 
R(xi,x) as long as o'L2\Xui^^Q. Dust will have little ef- 




Figure 8. Dust-free HI Frequency Redistribution Top: 
Reflection-line profiles for Xi = 5, 10, 20 are shown. The solid 
lines are simulations, while the dashed lines are the analytic pre- 
dictions -Ranlyt(2:), eqn. I37i . There are significant discrepancies, 
as discussed in the text. Middle: The reflection line profile for 
Xi = 10 is shown. The jagged solid line is the simulation, while 
the smooth solid line is the fit R{xi , x; a) , eqn. I39i , with o = 45/2 
and Xi given by eqn. 1381 . This simple rescaling of the analytic 
redistribution function gives accurate fits to the simulations. Bot- 
tom: The same as the middle panel, but shown on a logarithmic 
scale and over a larger frequency range. In all three plots, the 
sharp peaks at a; ±3 are an artifact of the acceleration scheme 
that we use; we do not calculate any core scatters, and photons 
which escape the core are all placed at x = ±3.3. Exact simu- 
lations indicate that there is indeed a "pile up" of photons just 
outside the core, but the peaks are not as sharp. In any event, 
the number of photons in the peaks are relatively insignificant. 



feet on frequency redistribution, except for extremely dusty 
or highly ionized clouds. 

For photons incident at frequency Xi on an optically 
thick slab, aro ^ 10"^, an analytic solution for the transmitted 
and re flected emission profile has b een derived bv_Neufeldl 
il990l) . extending the earlier work of lHarringtorJ Jl973l'l who 
obtained these results for the c ase Xi = 0. By taking the 
To ~* oo limit of eqn. (2.33) in iNeufeldl Jl990l) . we derive 
the analytic result for the reflected spectrum (normalized to 
unity) : 



-^anlyt(^l, ^) 



27r2 (a;3 „ 2.^3)2/6 +2.^4 



(37) 



When compared to simulations, -Raniyt(2;) is inaccurate in 
two respects. First, the actual peaks are shifted by an 
amount —2/xi (towards line center) compared to 7?aniyt- 
This can be compensated for by using a shifted incident 
frequency Xi in place of Xi, where 



2/x, 



(38) 
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Figure 9. Dusty HI Frequency Redistribution Top: The 
simulated reflection line profiles for Xi = 10 and three different 
t7°2i a-re shown: the thin line is ct^ji = Oi the light thick line 
is cr'^21 ~ 1 ^^'^ dark thick line is <t" 21 = 10- Bottom: The 
same as the top panel but on a logarithmic scale and over a larger 
frequency range. Over the frequency range < x < 15, the fre- 
quency redistribution function R{x, Xi) differs by < 20% between 
< fjiji < 1; hence, to lowest order it is independent of dust 
content for the cases we are interested in. 



Second, the peaks are significantly flatter than those given 
by -Raniyt- It is not sur prising that -Ran ivt is inaccurate, since 
the analytic results of iNeufeldl il99Ci) only hold when the 
photons traverse a large line center optical depth, aro ^ 10'^. 
The bulk of photons that reflect off a surface do not traverse 
such a large optical distance, and so the analytic result may 
not be accurate, which seems to be the case. However, one 
can obtain an excellent approximation to the redistribution 
function by modifying eqn. 13711 to include an additional 
fitting parameter a (see Appendix IBl for details): 



R{xi,x;a) 



3^2 



(39) 



Note that R{xi, x; a) is guaranteed to be normalized for all 
a. For example, -Raniyt is recovered with a = 6 and Xi — Xi . 
We find that our simulations for pure HI are well fit by 
a = 45/2, with Xi given by eqn. 1381 . To generate ran- 
dom frequencies x that obey the probability distribution 
R{xi, x; a), one draws a random univariate u G [0,1], and 
sets 



dx' R{xi, x' ; a) = F{x) 



(40) 



The frequency x is then given by functional inversion, x — 
F~^{u). Carrying out these steps on eqn. 13911 . done in Ap- 
pendix^ we find that the exiting frequencies x which obey 
the probability distribution R{x; a) can be generated by the 
equation 



X — [xi — xf\/a tan (ttu)] 
where u £ [0, 1] is a random univariate. 



(41) 



When dust is included, the profile peaks become slightly 
sharper and the tails fall off slightly faster. However, as 
shown in Figure |3] the pure HI distribution closely matches 
the dusty distribution when cTt2i < 10. In practice, we shall 
therefore always adopt eqn. 1391 with a = 45/2 for frequency 
redistribution in the line wing, since it is accurate except for 
galaxies with highly supersolar metallicities, which is un- 
likely at high redshift. 

When the incident frequency is in the Doppler core, the 
analytic fit, eqn. 1391 . breaks down, and the emission pro- 
file takes on an entirely different form. The emission profile 
roughly breaks down into two principal components: pho- 
tons that escape after only a few scatters, and photons that 
scatter enough times that they reach line center before es- 
cape. The former photons retain some "memory" of their 
incident frequency Xi, and produce emission peaks a.t x = Xi 
and X = —Xi. (the peak at —Xi is from photons that have 
undergone a single particle back-scattering off atoms with 
velocity V = XiV^""^ along the photon propagation direc- 
tion). The latter photons lose all memory of their initial 
frequency, and produce a broader emission peak centered 
on a; = 0. Accordingly, we fit the core redistribution func- 
tion R'^°'^°{xi,x) with three gaussians, centered on —Xi, 0, 
and Xi, respectively. Let us define the gaussian distribution 
with r.m.s. frequency a: 



G{x,a) 



2tv 



(42) 



By comparing to exact simulations, shown in Figure [TUl we 
find that a decent fit is given by 



R'^'^^ix^x) = 

[1-P{xi)] [lG{x + x,,A) + l 
+ Pixi)Gix,C) , 



G(a 



,B) 



where 

A 
B 
C 



0.7e" 
0.4 
0.5 
1.25 



(43) 



(44) 



The above fit works well when dust is included, since the 
effect of dust on photons in the Doppler core is small in 
absolute terms (see H3.1.11 . 

In summary, a simple analytic prescription for the sur- 
face frequency redistribution function R(xi,x) is to use 
equation 1391 with a = 22.5 and Xi = Xi — 2/xi when the 
incident photon is in the line wing \xi\ ^ 3, while equation 
1431 can be used when the incident photon is in the Doppler 
core \xi\ < 3. 



3.4 Surface Kinematics 

In this section we discuss two kinematic effects of surface 
scattering. First, we calculate the average scattering angle 
cosine g — (cos 6 scat), where Oscax is the angle between the in- 
cident and exiting direction. Second, we calculate the net fre- 
quency shift AV due to the Doppler shift induced by scatter- 
ing off a moving surface. In each case, we consider isotropic 
and perpendicularly incident photons, and use the surface 
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Figure 10. Dust-free Core Redistribution The simulated 
emission line profiles are shown for three different incident fre- 
quencies that lie within the Doppler core. As indicated by the ar- 
rows, the panels from top to bottome correspond to Xi = 2, 1, 0. 
In each panel, the thick line is the exact simulation data and the 
thin line is the fitting formula eqn. I43i . 



scattering angular distribution Dss{9), eqn. 12811 . The case 
of perpendicularly incident photons is of interest because it 
applies to photons bouncing around inside a spherical shell; 
most photons that strike a point on the surface last reflected 
off the far side of the shell, and hence incident photons have 
a strong bias to lay along the surface's perpendicular. 

We use the following conventions: the outward normal 
to the surface defines the z direction, an incident photon 
has direction ni such that fia ■ h\ < with polar angle 

61 and azimuthal angle 01, an exiting photon has direction 
n2 such that ns • ^2 > with polar and azimuthal angles 

62 and (f>2. For example, in Cartesian co-ordinates, fii = 
(cos sin ^1 , sin 01 sin ^1 , cos^i). The surface has a bulk 
velocity Vs, with a perpendicular component V±_ = fia ■ Va- 
For an isotropic angular distribution of incident photons, the 
polar angle distribution of incident photons is DisoiOi) = 
sin 01, which is normalized to unity over di G [7r/2, vr], and 
4>2 is uniformly distributed over 2ii. For exiting photons, 
eqn. H28^ . the polar angle distribution is given by Dss{02) = 
2 sin 02 cos O2 , which is normalized to unity over 62 G [0 , 7r/2] , 
and 02 is uniformly distributed over 2n. 



the angles is 

5iso = 



/'27r /'7r/2 /'27r 

d^i / d0l / d6l2 / d02 X 

tt/2 Jo Jo Jo 



(45) 



\^D,so{9i)Dssie2) hi ■ h2^ 



Azimuthal symmetry eliminates all the terms in h\ ■ 112 
that have an x and y contribution, leaving just the term 
cos 611 cos 62. The integral over 9i and 62 separate, giving 



(46) 



The same steps can be carried out for perpendicularly inci- 
dent photons, ni = —z, resulting in 

5x = -|. (47) 

Surface scattering is characterized by a net average back- 
scatter, 3 < 0. 



3.4-2 Frequency Shift Due to a Bulk Surface Velocity 

If the surface has a bulk velocity, then the frequency of the 
scattered photons suffer a net Doppler shift, due purely to 
the surface motion. Consider a photon with frequency V 
striking a moving surface. In the surface rest frame, the 
photon has incident frequency V' = V — n\ ■ Vs- An ex- 
iting photon with frequency V" in the surface rest frame 
has an exiting frequency V'" = V" + h2 ■ Vs in the original 
("lab") frame. Therefore the surface motion induces a net 
frequency shift 



AKm = [hi 



n2 



Vs 



(48) 



which is in addition to any frequency shift from scattering 
within the surface. Averaging over an isotropic incident an- 
gle gives 

1 . 



{hi)iso 



(49) 



For perpendicularly incident photons, we simply have (ni) — 
—z. For the exiting direction, averaging over Dss{d2) gives 



{h2) = -z 



(50) 



Thus, the average frequency shift for isotropic and perpen- 
dicularly incident photons are, respectively, 



7 

{AV;m)iso = -v± 
6 



(51) 
(52) 



where, as stated above, we define V± = hs ■ Vs = z ■ V3. li 
the surface is moving away from (towards) the incident pho- 
tons, V± < 0, then surface scattering causes a net redshift 
(blueshift) of the photon frequency. 



3.4-1 The Average Scattering Angle 

The scattering angle cosine, also called the scattering asym- 
metry parameter, is defined hy g = (cos^scat) = {hi ■ 712). 
For an isotropic incident angle distribution, the average over 



4 ANALYTIC MULTI-PHASE TRANSFER 

In this section, we build an analytic model for estimat- 
ing the escape fraction and line width for Lya escaping 
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from a multi-phase region composed of dusty, optically thick 
clumps. Compared against Monte-Carlo simulations, these 
analytic estimates give remarkably accurate results. We con- 
sider both stationary clumps and clumps with a Maxwellian 
bulk velocity distribution, but postpone the discussion of 
bulk gas outflow/inflow until H5.2I We derive a general pho- 
ton escape fraction formula in terms of two parameters: the 
mean number of surface scatters in the absence of absorp- 
tion Ao, and the cloud albedo Ec, for the case of coherent 
scattering. The parameter Ao is purely geometric and in- 
dependent of photon frequency, as long as clouds are very 
optically thick. We derive fits for A/o for a variety of multi- 
phase geometries of astrophysical interest, and compare the 
resulting analytic escape fractions to Monte-Carlo simula- 
tions of coherent scattering. 

Applying the coherent scattering formulas to Lya radia- 
tive transfer requires some care, since in this case is fre- 
quency dependent through the frequency dependence of the 
Lya scattering cross-section. The cloud albedo, therefore, 
changes as the photon executes a random walk in frequency. 
Nonetheless, we find that the analytic formula for coherent 
scattering can be extended to Lya scattering as long as tc is 
evaluated at the characteristic escape frequency. In H4.3I we 
study how the characteristic escape frequency depends upon 
the frequency redistribution per surface scatter and the ran- 
dom bulk motion of the clumps. We derive estimates of the 
escaping line profile r.m.s. width and FWHM, and compare 
these to simulations of repeated surface scatterings. 

4.1 Escape Fractions 

A more detailed prediction for the escape fraction can be 
given than the simple scaling laws in ^ We derive a generic 
escape fraction formula in terms of the probability of absorp- 
tion per scattering, e, and the average number of scatters for 
escape in the absence of absorption, Ao. In the context of 
multi-phase transfer, each "scatter" refers to an entire sur- 
face scattering, in which case e is the cloud albedo: e = ec- 
Central to the analysis is the average number of interac- 
tions for escaping photons. A/". Let e be the average absorp- 
tion probability per interaction, and define D{n) to be the 
probability distribution for photons to interact n times be- 
fore escaping, when e = (i.e., in the absence of absorption). 
For a constant e, the escape fraction can be written 



From eqn. 155L A/fj = A/'(0) is given by 



/e = ^(l-6)"i3(. 



(53) 



The average number of interaction for escaping photons can 
likewise be written as 



^ oo 

AA=-^(l-.)"n75(n) . 

n—O 

From these two expressions, we derive 



(54) 



Ai) = lim — — In /e 

e^o de 



(57) 



To derive escape fractions in an arbitrary geometry, let 
us consider two limits: narrowly beamed forward scatter- 
ing, and front-back symmetric scattering. The latter refers 
to the g = (cos^sct) ~ case, where there is equal prob- 
ability to scatter in the forward and backward directions; 
both isotropic and dipole scattering satisfy this, though dust 
grains (gd = 0.5) and cloud surfaces (piao = —1/3, from eqn. 
1461 1). do not. As we shall show, however, the escape fraction 
for scattering for any value of g < 1 can be based upon the 
g = case as long as there are sufficient scatters that the 
photon's trajectory is randomized. 

The case of purely forward scattering is trivial: since 
the photon does not change direction, we have 



(58) 



where we have used eqn. 15711 in the final step. Note that the 
escape fraction is same as if there were no scatterers, since 
scattering does not alter the photon trajectory. For front- 
back symmetric scattering {g — 0), the escape fraction is 
(see Figure mn in Appendix IHt: 



cosh(ei/Vr2 + 2r) cosh (V2&A/^) 



(59) 



where we applied eqn. 1571 to obtain A/t) ~ \ (t^ + 2t) in the 
second step^. Although this formula for the escape fraction 
is derived assuming gr = 0, we show in Appendix |0 that 
this formula is valid for any g as long as there are enough 
scatters to randomize the photon's direction. Specifically, 
if No > ri'{g) (where n'{g) is given by eqn. ijCS^ V then 
eqn. 15911 can still be used. Otherwise, eqn. 15811 is a bet- 
ter approximation to the escape fraction — since a situation 
with such few scatters and/or a strongly peaked forward 
scattering profile converges to the straight-line trajectory 
case. Note that in both these limits, the escape fraction 
only depends upon a single parameter, eAo. The average 
number of interactions required for a photon to be absorbed 
is Ma ~ so the controlling parameter is equivalent to 
eAAo =A'o/A'„. 

The average number of interactions for an escaping pho- 
ton, A/", can be derived by applying eqn. 1551 to the appro- 
priate escape fraction formulae, eqn. 15811 or eqn. 1591 . For 
straight-line trajectories, the result is 



A" = (1 - 6)A'o , 

while for the random walk trajectories 
tanh {V^dT^) 



A" = (1 - e)A'o ■ 



A/2iAl^ 



(60) 



(61) 



-(1-^)3- In ./c . 
de 

This can be inverted to express /c in terms of A/": 
/.=exp - / d.^ 



(55) 



(56) 



Within the Eddington approximation, the escape fraction from 
a source in the mid-plane of a slab is /o = 1/ cosh ^v^S e^^^rj, 
as derived, for example, in iNeufeldl ^^^). Instead, we propose 
eqn. I59i . which from our simulations is more accurate for 1-D 
scattering (see Appendix Ict. 
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In the limit v^eA/^ < 1, we find N ^ No - eA/'o(|A/'o + 1), 
and in the hmit y/TFJo 2> 1, we have TV -^^^ ^SeA/o- 

The entire effect of the cloud geometry is characterized 
by a single number, Mo - Directly calculating Ao from a given 
geometry is not practical except for the simplest geometries. 
In general, given a clump geometry, Ao must be computed 
via a simulation, where one can use the exiting angle distri- 
bution, eqn. 12811 . for each surface refiection. In practice, we 
find that for many generic geometries expected to crop up 
in astrophysical applications, an appropriate "line of sight" - 
averaged Ao can be accurately determined from simple ge- 
ometric parameters (such as the cloud covering factor /c). 
We now proceed to do so. 

4.2 Example Geometries 

We now test the accuracy of the analytic formula eqn. I59II . 
and investigate its application in a variety of multi-phase 
geometries. For each geometry, we fit Ao as a function of 
the appropriate natural geometric parameter (such as cloud 
covering factor fc) using Monte Carlo simulations. We then 
calculate the escape fraction using A/o and eqn. 1591 . and 
compare this to simulations. In several cases, we find that we 
obtain better fits by rescaling with an order unity fitting pa- 
rameter K, where we use kecAo in place of EcAo in the escape 
fraction formulae eqn. (1591 . In general, even with no correc- 
tion factor, eqn. I59II is accurate to <, 10% when /e ^ 10%, 
and is generally correct to within a factor 2 for /e < 10%. 
When the escape fraction is very small, the photons that 
do escape comprise the rare trajectories. As their transfer 
behavior can depend sensitively on the specifics of the ge- 
ometry, it is not surprising that eqn. 15911 breaks down when 
/e is very small. In any case, Lya emission is undetectable 
for such small /c, so these cases are of little observational 
consequence. 

Some notes about our Monte-Carlo simulations; for sim- 
plicity, we assume that the region in between the scattering 
surfaces is empty (we justify this assumption in ii4.5ll . We 
also assume that scattering surfaces are extremely optically 
thick, so that the surface scattering approximations of [0 
apply. In particular, when a photon hits a gas surface, it has 
a net probability ec of being absorbed; if it survives, then its 
exiting angle relative to the surface normal follows the dis- 
tribution Dss{0) = sm20, and the exit location is the same 
as the point of incidence. 

4-2.1 Spherical Clumps 

The canonical example of a multi-phase geometry is a spher- 
ical region populated by randomly placed, optically thick 
spherical clumps. Such clumps tend to be cool-phase gas 
such as molecular clouds, which arose via thermal instabil- 
ity. The natural geometric parameter is the mean number 
of clouds intersected along a random line of sight, called the 
cloud covering factor fc- The covering factor is analogous 
to the interaction optical depth r for homogenous media. 
For a central source, fc is measured from the region center 
to the edge. We computed A/o for various covering factors 
and clump radii distributions. As shown by the solid line in 



the top panel of Figure ITTI we find that a general fit for any 
clump radii distribution is 

AAo = /c' + pc - (62) 

This highlights the fact that when surface scattering applies, 
the spherical clump geometry does not depend upon the size 
distribution of the clumps nor their volume filling fraction; 
the entire geometry is chara cterized by a si ngle parameter, 
fc- This was postulated bv 'Ne ufeldl l)l99lfl : we have con- 
firmed this insight in our simulations. The scalings in eqn. 
1621 1 can be compared against the usual random walk for- 
mulae, A/" ~ r (for T ^ 1), and JV ~ (for r ^ 1), for 
scattering with front-back symmetry g = (e.g., Rybicki 
& Lightman 1979, p. 35); here, fc plays the role of r. Our 
scalings are slightly different, since g = —1/3 for our clouds. 
For 1-D scattering, Ao ~ 5 (t^ + 2r) (shown as a dotted line 
in Fig, mil . By comparison against eqn. 1621 . the spherical 
clump model is analogous to 1-D radiative transfer with the 
substitution r ~ V^fc - In the bottom panel of Figure 1111 
we show that the escape fraction formula, eqn. I|59^ . works 
well for a variety of covering factors for constant ec- 

While the covering fraction is an unknown free param- 
eter, it is easy to see why /c ~ 1 is reasonable. Suppose 
the cold clouds constitute a mass fraction /m of the galaxy, 
and are overdense by a factor 6 relative to the intercloud 
medium. Assuming pressure balance, 5 ~ Ticm/Tc ^ 100, 
for Tc ~ lO^'K and Ticm ~ lO^K for the cloud and inter- 
cloud medium temperatures respectively. The volume filling 
factor of clouds is then fv = fhi/S- The number density 
of clouds is ric ~ /v/K:, where Vc ~ is the volume of a 
typical cloud. The cloud covering factor is: 



4-2-2 Random Surfaces 

An abstraction of the spherical clump geometry is to have 
the photons strike a surface after traveling a distance I, 
where I is drawn from a given probability distribution, and 
where each surface is randomly oriented with respect to the 
photon direction. We have investigated the exponential dis- 
tribution exp{— fc£/R), where R is the radius of the region 
and fc is the covering factor. The photon escapes once it 
leaves the spherical region of radius R- When a photon trav- 
eling in the direction hp strikes a surface, the outward nor- 
mal of the surface hs is drawn from an isotropic distribution 
such that hp ■ hs < 0. As shown by the top panel in Figure 
1121 a good fit for the scattering number is 

Afo = hc^ + fc (64) 
5 

The random surfaces model is faster to simulate, and any 
results reliably apply to the spherical clump model when 
expressed in terms of Ao . The bottom panel of Figure 1121 
compares the escape fraction formula eqn. 1591 to simula- 
tions. 

This procedure for generating random surfaces can in 
fact be used to quickly simulate any arbitrary geometry, with 
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Figure 11. Random Spheres Geometry Spherical clouds of 
radius r randomly populate a spherical region with radius R. 
The circles have clouds with the same radius, r/R = 1/20. Top: 
The number of scatterings A/q as a function of the cloud covering 
factor fc- The squares have clouds with the same radius r/R = 
1/40. The diamonds have a distribution of cloud sizes with the 
distribution D{r) oc 1/r^, with a minimum cloud radius r/R = 
3/200 and a maximum cloud radius r/R = 1/20. The solid line 
is the fit eqn. 1621 . while the dotted line is the naive 1-D relation 
J\fQ = i(/^2^_2j^)^ extrapolated from A/q = ^{t^+2t), and t 
fc (in fact, t —> \/2/c is the appropriate substitution). Bottom: 
Photon escape fraction /e as a function of cloud albedo ec- From 
top to bottom, the covering factors are fc = 0.48, 1.34, 2.14, 3.0. 
The lines are based on eqn. 1621 . 



a suitable characterization of the probability distribution of 
path lengths £ and surface orientations ha, though, of course, 
the fitting formulae for Ao will depend on these quantities. 

4.2.3 Shell with Holes 

Another geometry that frequently arises in astrophysics is 
that of a shell of material surrounding a photon source — 
for instance, in stellar and galactic outflows. Much work 
has be en done on Lyg rad i ative t ransfer through opaque 
shell s llTeno rio-Tagle et al I Il999l: lAhn. Lee, fc Led 120031: 
I Ahn 111.2004) . but the effect of gaps in the shell has not been 
investigated. Since a completely homogeneous shell of gas 
is rarely, if ever, observed, a shell with holes is an interest- 
ing geometry. The natural geometric parameter here is the 
fraction of the solid angle covered by the shell, fcov (i.e., the 
gaps comprise a total solid angle 4n{l — /cov)). To estimate 
Ao, assume that during each bounce the photon has a prob- 
ability of (1 — /cov) to escape through a gap. This leads to 
the expression 

Afo ~ /cov /(I - /cov) , (65) 




Figure 12. Random Surfaces Geometry Top: A'o vs. fc. 
Points arc simulation results, while the solid line is the fit eqn. 
I64i . As in the top panel of Figurc lTTl the dotted line is the naive 
1-D power law A/q = ^{fc^ + '^fc)- Bottom: /c vs. ec. From top 
to bottom, the covering factors are fc = 1, 3, 5, 10. The solid 
lines use EcA/q in the escape fraction formula eqn. I59i . while the 
dashed line uses fcecA/o, with a correction factor re = 1.29. 



which is shown by the solid line in the top panel of Figure 
1131 We computed A/q for shells with a random placement of 
non-overlapping circular gaps, and found that many small 
gaps were equivalent to fewer large gaps with the same /cov, 
as to be expected. As shown in the bottom panel of Figure 
1131 the analytic form for /e does well where expected, and 
only breaks down when /cov is near unity and ^ 0.3. This 
agreement is quite remarkable, given that the escape fraction 
formula (I59II was derived assuming a very different scatter- 
ing particle geometry (homogeneous media). This gives us 
confidence that the impact of geometry on the escape frac- 
tion can indeed be encapsulated by a single parameter ecAo. 



4.2.4 Open Ended Tube 

To illustrate how well photons can escape through cracks 
and gaps in optically thick material, we consider photons 
escaping from the middle of an open-ended tube. This geom- 
etry may apply to star forming regions or AGNs where out- 
flows have punched (bipolar) cavities through a surrounding 
gas cloudy allowing the photons t o escape thr ough the cavi- 
ties (e.g. jShoDbcU & Bland-Hawt horn I lll99a) ). although we 
do not investigate the dependence on the opening angle. The 
natural geometric parameter is the ratio L/R, where L is the 
total length of the tube and R is the tube's radius. The av- 
erage length traveled along the tube per scatter is £ ^ R. 
Therefore A/o is estimated from L/2 \/ Mo£, which implies 
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Figure 13. Shell with Holes Geometry Top: The scattering 
number Afo as a function of the sheU covering factor fcov ■ Simula- 
tions were run for shells with one gap (light circles) and five gaps 
(dark circles), which fall almost exactly on the fit eqn. 1651 . As 
expected, the scattering number Afo only depends on fcov, and 
not the number of gaps. Bottom: /o vs. tc- From top to bottom, 
the covering factors are fcov = 0.5, 0.7, 0.8, 0.9. The lines show 
the fitting formula eqn. I59i . with kecA/q as the variable, with 
K = 1 (solid lines), and k = 0.8 (dashed lines). 



Ao ~ (L/R)^. As shown by the solid line in the top panel of 
Figure IT^ a decent fit for the number of surface scatters is 

J^o = ^iL/Rf + ^{L/R) (66) 

As shown in the bottom panel of Figure [T^ the escape frac- 
tions fit the data very well, except as — > 1- As previously 
discussed, for such low escape fractions, rare escaping pho- 
tons follow unusual trajectories not well captured by our 
formalism. In any case, for such low escape fractions, Lya 
emission cannot be observed, and the results have no obser- 
vational relevance. 

4.3 Line Widths 

In this section we shall consider two separate effects on the 
Lya frequency as it escapes: the frequency redistribution 
due to thermal motion of atoms, as well as random bulk 
motion of the scattering surfaces. The effect of global out- 
flow/inflow on the line profile is discussed separately, in H5.2I 
We consider two simple diagnostics of the profile: the r.m.s. 
frequency shift E and the FWHM F. We ran simulations 
of dust-free surface scattering to determine E and F as a 
function of the number of surface scatterings, Uss- Although 
escaping photons vary in the number of scattering surfaces 
they encounter before escape, in practice we find that line 
profiles can be accurately characterized by the average num- 




Length/Radius 




Figure 14. Tube Geometry Top: The solid line is the fit for 
A/q, eqn. 1661 . and the circles are simulation results. Bottom: From 
top to bottom, the length to radius ratios are L/iJ = 2, 5, 10, 20. 
The solid lines are given by eqn. I59i . with no correction factor 

iK=l). 



her of scattering surfaces encountered. We now derive formu- 
las for E(nss) and T{nss) for resonant scattering frequency 
redistribution and random bulk surface velocities. 



4.3.1 Resonant Scattering 

By using the analytic approximation to the surface fre- 
quency redistribution function R{xi,x) derived in H3.3.5I 
we simulated repeated surface scattering for an ensemble 
of photons. Carrying this out using the exact Monte-Carlo 
simulation is not computationally feasible. For a few cases 
we compared the results of the analytic approximation with 
an exact Monte-Carlo simulation and a Monte-Carlo simu- 
lation that approximates the core scatterings ( H3.1.1II . and 
all three are in good accord. This gives us some confidence 
that the inaccuracies in the analytic approximation do not 
compound significantly for the number of scattering surfaces 
we investigated. As shown in Figure IT51 decent fits for the 
resonant scattering E and F are 

E.s(n,,) = ^nss V'"' (67) 

1 + Uss 

The behavior can be easily understood. If atomic thermal 
motions are responsible for frequency redistribution, then 
the line profile quickly relaxes to a Gaussian whose FWHM 
is determined by the width of the Doppler core. However, the 
profile also acquires non-Gaussian tails from rare scattering 
events; these tails dominate the r.m.s. line width E, hence 
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Figure 15. Line Widths Top Left: The r.m.s frequency shift 
S (filled circles) and the FWHM F (open circles) are shown as 
a function of the n**" surface scatter Uss for pure atomic scat- 
tering. The lines are eqn. 1671 . Top Right: The r.m.s frequency 
shift E and the FWHM F are shown as a function of the n"" 
surface scatter for surface motion with a Maxwellian distribution 
with r.m.s. speed V^. Bottom Left: The r.m.s frequency shift S 
as a function of the n"" surface scatter when both atomic ther- 
mal scattering and random cloud motion effects are combined. 
The four different symbols correspond to four different values of 
yc ^ydop. ^j^g tj-jangles, circles, squares, and diamonds correspond 
to y<=/\/dop ^ ;^o, 5, 3, 1 respectively The lines are eqn. IS^ . 
Bottom Right: The FWHM F as a function of n^s; otherwise all 
symbols retain the same meaning as the bottom left panel. 



E oc Uss- For resonant scattering only, the FWHM F is a 
more accurate measure of the line profile than E. Examples 
of the spectra after repeated surface scatters is shown in 
Figure [TC] 



4.3.2 Surface Motion 

As shown in H3.4I when photons scatter off a moving sur- 
face there is a net frequency shift per surface scatter of 
{AV} ~ V±, where V± is the velocity along the outward 
normal. If the moving surfaces have a Maxwellian velocity 
distribution with r.m.s. speed V^, then we expect that the 
induced r.m.s. frequency shift after Hss scattering surfaces 
will scale as Esin(fiss) ~ y'riss V^. For a Gaussian distri- 
bution, the FWHM and the standard deviation are related 
by F « 2.2E, and so we expect Fsni(w3s) ~ 2.2y/n^ V^. To 
check this, we simulated repeated surface scattering assum- 
ing an isotropic incident angle distribution and the surface 
scattering exiting angle distribution, eqn. I|28|l . We indeed 




Figure 16. Repeated Surface Scatters: Atomic Motion 
Only The normalized frequency distribution after repeated 
scatters off a stationary, dust-free, surface. Shown are nsa = 
1, 5, 10, 15, which correspond to increasing widths in the plot. 
An ensemble of photons are initially injected at line center, and 
their frequencies are tracked after repeated cloud scatters. For 
each cloud scatter, the analytic formula of H3.3.5l is used to gen- 
erate the photon's exiting frequency. 

find that 

Esm(?iss) = V (68) 

Fsm — 2.2^ Tiss V , 

which is shown in the middle panel of Figure IT^ 
4.3.3 Combined Effect 

In the two bottom panels of Figure 1151 we show the com- 
bined effects of resonant line broadening and surface motion. 
For most multi- phase geometries, 3> V^"^, so line broad- 
ening is dominated by cloud motions. In this regime, the 
line width can be accurately estimated by eqn. 1681 . Only if 
cloud motions are small and comparable to atomic thermal 
motions, does the behavior change. In this case, 

the FWHM does not increase with the number of scatters; 
instead it "thermalizes" to the characteristic Doppler width. 
The line profile is more accurately described by eqn. 16711 . 
When V''>V'^°'^, the total r.m.s. width is (note that the 
dispersions add linearly rather than in quadrature) 

E(nss) = 2Ers(ni,s) + Esm(".ss) ■ (69) 

The FWHM has a more complex behavior because the 
Doppler core tends to retard the increase in F beyond the 
size of the Doppler core. 

4.4 Lya Escape Fractions 

In this section we show how Lya multiphase transfer can 
be handled analytically based on the results in the previ- 
ous subsections. We shall first estimate the typical escape 
frequency of Lya photons Xc, which provides the typical 
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surface absorption albedo ec(xe)- With this estimate of ec 
and Ao, we can calculate the typical Lya escape fraction /e 
using eqn. H59^ . 

We find that an adequate approximation for the typi- 
cal escape frequency is given by a slight modification of the 
r.m.s. line width E(nss), eqn. (Iti91 . where the surface mo- 
tion can be either a random Maxwellian cloud velocity with 
r.m.s. speed or a bulk outflow with typical outflow speed 
(see H5.2I below). It is not correct to approximate the 
typical number of scattering surfaces Uas with the average 
number of scattering surfaces in the absence of absorption, 
Mo, since in general n^s ^ Ao once absorption is taken into 
account. A more appropriate measure of Uss is given by N , 
eqn. II61II . the average number of scattering surfaces encoun- 
tered before escape for a fixed cloud albedo Ec- However, 
since ec is frequency-dependent, we still need to estimate 
the escape frequency x^. We do so iteratively. We first esti- 
mate Xa assuming that absorption is unimportant, 

xl = E(A'o) , (70) 

which obviously overestimates Xc- We can then estimate M 



tanh V2e?Ao 



where 



(71) 



(72) 



with E(nss) given by eqn. H69|l . 

With the above estimate of Uss ~ A/", a better approxi- 
mation to Xc is given by the next iteration. 



x', = E(AA) 



(73) 



which gives a better approximation for the typical cloud 
albedo 



1 + 



(74) 



At this point one could iterate again — using e'^ to find a 
better approximation to riss — but we find that stopping 
after one iteration provides escape fractions in good accord 
with simulations. From eqn. 159L the escape fraction given 
by 



(75) 



cosh x/2ecAo 

In Figure 1171 we compare this analytic escape fraction 
to simulations of radiative transfer through the random 
surfaces geometry (see H4.2.21 for both random cloud mo- 
tions and a bulk cloud outflow, and for several amounts of 
dust. The bulk cloud outflow is purely radial, with the same 
speed at all radii, which approximates galactic winds out- 
side the initial acceleration zone. As can be seen, the an- 
alytic approximation captures the simulated escaped frac- 
tion to ~ 20% when cloud motion dominates atomic ther- 
mal motion, V^'^ > 3V''^°p ~ 40km/s. As explained in 
once the effects of the Doppler core become important, the 
r.m.s. frequency shift E is no longer a good measure of the 
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Figure 17. Analytic Lya Escape Fraction Examples The 

Lya escape fraction fc as a function of the typical speed of the 
bulk cloud motion, V^. Simulations of Lya transfer through the 
random surfaces geometry I' il4.2.2l with fc=3 were run for two 
types of bulk cloud motion: a Maxwellian velocity distribution 
(open symbols), and a purely radial outflow (filled symbols). In 
the former case the x-axis corresponds to the r.m.s. speed of the 
clouds, while in the latter case the x-axis corresponds to the out- 
flow speed, which is the same at all radii. Two different dust 
contents were simulated, both with e^j = 0.5: = 2 (circles) 

and o'!^2i ~ 0.2 (squares). In all cases the gas temperature is 
10^ K. The analytic approximations (lines) are generated by the 
steps outlined in text. 



typical frequency: E will overestimate the typical escape fre- 
quency, and hence lead to an overestimate of the absorption. 
This is seen in the flgure for <, 40km/s. However, since 
the amount of absorption is typically not significant when 
the escape frequency is 40 km/s, for most purposes using 
E allows one to accurately estimate f^. Note that for bulk 
outfiows with the same characteristic speed V^, the escape 
fraction is smaller than for random cloud motion. The line 
profile for random cloud motion is approximately Gaussian 
centered on the line center, with a standard deviation of V. 
In contrast, the line profile for a bulk outflow is has a mean 
at V^. For a given V^, a bulk outflow produces more photons 
far from line center than random cloud motion, and hence 
the bulk outflow causes more Lya absorption. 

In summary, the escape fraction depends upon flve pa- 
rameters: the typical cloud speed V, the number of surface 
scatters in the absence of absorption A/q, the gas tempera- 
ture T, and the dust parameters ct'* and e''. 



4.5 Dust and Gas Between the Clumps 

Thus far, we have only considered radiative transfer off 
opaque surfaces, and ignored absorption and/or scattering 
in the optically thin hot intercloud medium (ICM) . We treat 
the hot intercloud medium as having a very low neutral 
HI fraction (for gas in coronal equilibrium at T ~ lO'^K, 
xhi ~ 10~^), as well as being dust-depleted, due to sputter- 
ing and other dust-destruction processes. ICM resonant scat- 
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tering and absorption is easily incorporated in our Monte- 
Carlo simulations at little computational cost, and we have 
experimented with various prescriptions for the ICM in our 
runs. In most cases, we find that radiative transfer within 
the ICM can be neglected, and here we show some simple 
estimates demonstrating why this is the case. 

Let us first consider dust absorption in the ICM, as- 
suming that resonant scattering in the ICM is negligible. It 
is easy to show that for a multiphase medium, the relative 
fractional column densities of a species i in phases X, Y is 
simply given by the relative mass fraction of species i in that 
phase: 



NL 



fli 



fM,X 
fM,Y 



(76) 



where x^x is the mass fraction of species i in phase X. Thus, 
for instance, N^i/N^"^ ^ fMA^m"^]'^ » 1 for x^g^ < 1, 
and the observed HI column density N^\i ~ 1 is strongly 
dominated by gas in the cold phase. 

We can use this to estimate scattering in the ICM. Due 
to its random walk as it scatters off optically thick surfaces, 
a Lya photon traverses an ICM column density ~ vCV times 
larger than for a straight line path. At a characteristic escape 
frequency Vy ~ 1, it therefore encounters an ICM HI optical 
depth: 



Thi = 0.07iVHl,2 



/m.icm 
0.7 



5 



ICM 
10-4 



Hence, resonant scattering in the ICM is negligible. It only 
becomes important when the photon is within the Doppler 
core V2 ~ 0.26 for the parameters above). However, even 
in this rare case ~few, comparable to the number of 

surface scatterings A/". 

What about dust absorption in the ICM? Let us sup- 
pose that due to dust depletion, only /^"""^ ~ 0.05 of all the 
dust in the galaxy is in the ICM. If the total dust optical 
depth through the cold phase is ~ 1 A^hi'2i''"-21, then the 



ICM 



Id 



ICM„ 



total dust optical depth through the ICM is 
The mean free path of a Lya photon is I ~ rgai / VJV. Hence, 
between each bounce, a Lya photon traverses an optical 
depth r^°""<:= ~ rf^'^/y/M, and has a probability of ab- 
sorption in the ICM of Pl^^^^ « 1 - e 
or: 



_bouncc 



0.02 



'1/2 



Id 



ICM 



0.05 



(78) 



By contrast, during each surface scatter, the photon has an 
absorption probability of: 



,1/2 



(79) 



Thus, photons are more likely to be absorbed on cloud sur- 
faces, rather than in the ICM, justifying our neglect of ICM 
dust absorption. 

Obviously, all of these statements are parameter depen- 
dent and there are cases when scattering and absorption in 
the ICM cannot be neglected (if the IGM dust or HI con- 
tent is high). For instance, if the HI fraction in the ICM 
is high, then Lya can be strongly quenched. This may be 



partly responsible for the variation in the Lya equivalent 
widths amongst different galaxies (see discussion in >I5.H . 



4.6 Accelerated Radiative Transfer on a Grid 

The approach of this paper is to identify all optically 
thick surfaces in a Monte-Carlo simulation and apply the 
scattering properties identified in |3]to them, thus afford- 
ing vast computational speed-ups. This should be readily 
amenable to radiative transfer in numerical simulations. For 
this scheme to be accurate: (i) the surfaces must be suffi- 
ciently optically thick that transmission is negligible; i.e., 
they should satisfy equation 1321 . (ii) The approximation 
that the photon is either reflected or absorbed on the spot 
without significant spatial diffusion must hold. We now dis- 
cuss this second requirement. 

For the "on the spot" approximation to hold, the pho- 
ton's mean free path li should be significantly smaller than 
the grid cell size Lcdi. The photon typically moves a distance 
~ y/Jf^ii ~ Ui (see ja.l4l for discussion of AA="'f ) whilst 
scattering within the surface. Thus, we require liija Lccii 
where a ~ 1/10 is a constant that designates the desired 
level of accuracy for the approximation. The "on the spot 
approximation" is accurate if the HI density is larger than 
a critical density: 



HI 

'^ccll 



>n = 



13 



(a/0.1) (Lccii/pc) 



(80) 



where is the incident Lya frequency shift off line center 
in the rest frame of the HI in the cell (in velocity units). 

Alternatively, if one is willing to sacrifice some spatial 
resolution, one can consider a group of neighboring cells 
which are opaque across their total length (and thus sat- 
isfy equation 1)32^ l even though an individual cell may not 
be opaque. For cubic blocks of A^bik cells per side, the entire 
block surface will act like an absorbing mirror if 



/ \HI 



\ * 

n 



13 



Vn 



(a/0.1) (A^bikLccii/pc) 



(81) 



where (n)bik is the average HI density within the block. The 
surface approximations break down if the block is strongly 
inhomogeneous, i.e., the mean free path changes over a 
length scale that is much shorter than the block length. 
This is equivalent to |nccii/V(ncoii)| ^ /JA'^bikicoii, where 
/3 ~ 1/10 is a constant that designates the desired level 
of accuracy for the approximation. In terms of cells on a 
grid, let n{2\^\\ and n(l)ceii be any two neighboring cells 
in the block. The second condition on the absorbing mirror 
approximation takes the form 



n(2)c 



n(l). 



- 1 



(82) 



We look forward to implementing these ideas in numer- 
ical simulations in the future. 
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5 APPLICATIONS 

We briefly discuss two examples of applications of our ra- 
diative transfer framework. Many more extensive studies are 
possible. 



5.1 Lya Equivalent Widths 

Most Lya photons are produced in the HII regions sur- 
rounding sources of ionizing radiation, where roughly 2/3 of 
the ionizing photons are converted into Lya photons (under 
case B recombination). Hence, in the absence of radiative 
transfer effects, the equivalent width measures the number 
of ionizing photons emitted relative to the UV continuum 
near 1200 A. There are numerous examples of high-redshift 
sources which have equivalent widths which are too large 
to be produced by conventional stellar populations. About 
~ 2/3 of the SCUBA submm galaxies with accurate posi- 
tions from radio detections have Lya in emission , many with 
equiv alent widths too great for stellar sources llSmail et alJ 
2004). The mvsterio us Lya emitters at 2: ~ 3.1 observed by 
ISteidel et afl (I2OO0I) have enormous Lya fluxes, but no ob- 
served continuum. Finally, the LALA survey detects many 
high redshift z = 4.5, 5.7 sources with equivalent widths 
EW ^ 150 A signi ficantly in excess o f any known nearby 
stellar population jRhoads et aljl2003f) . An AGN origin is 
unlikely because follow-up observations show no signs of 
the X-rays and high-ionization line s expected for a ty pe II 
quasar source dWang et al. 120041 : iDawson et al. l2004) . An- 
other possibility is that the Lya emission is due to an ex- 
tremely top-heavy population of massive PopIII stars. How- 
ever, there are no signs of the strong Hell omission at 1640 A 
expected from metal- free stars (JDawson et al. 2004). 

Anot her poss i bility for high Lya EWs, originally sug- 
gested bv lNeufeldl il991f) . is radiative transfer effects. If the 
continuum is more absorbed than Lya photons during the 
escape from the host galaxy, then the equivalent width of 
the transmitted spectra is larger than the equivalent width 
of the source. The initial and transmitted equivalent widths 
are basically related by the ratio of Lya to continuum escape 
fractions, 

EWout ~ ^ EWsrc , (83) 

where EWsrc is the source equivalent width and EWout is 
the equivalent width for the escaping photons. In order for a 
"normal" starburst IMF with an intrinsic equivalent width 
of EWarc ~ 150 A to produce an observed equivalent width 
of EWout ^ 300 A, then radiative transfer must account for 
a "boost" by a factor of at least 2 — 3. For sources where no 
continuum is observed, the continuum must be preferentially 
extinguished by an even larger factor. 

Let us now estimate equivalent width boosts in our 
multi-phase model to see if this is possible. For any multi- 
phase medium where the gas resides in clumps that are very 
opaque to Lya , the surface scattering approximations apply, 
and so the Lya escape fraction can be analytically estimated 
as in H4.4I What about the continuum escape fraction? For 
simplicity, assume that each gas clump is not opaque to dust 



extinction: <,1, where is the dust extinction (scatter- 
ing-!- absorption) optical depth across a clump diameter*. 
Since the self-shielding effect of clumpy gas is therefore small 
for the continuum photons, the effective dust distribution 
is approximately homogenous for the continuum radiative 
transfer. The escape fraction for a photons injected in the 
middle of a homogenous medium, with an absorption albedo 
e'* ~ 1/2, is approximately that given by eqn. 159|l . 

/f"^^l/cosh[y?^lM2 + 27^] , (84) 

where t'' = N a^^ is the average dust extinction optical depth 
through a region with average HI column density A'^. Figure 
ITsl shows how the ratio /c^^°//f ™ varies as a function of A^o 
for a fiducial set of multiphase gas parameters. A substantial 
"boost" in the transmitted equivalent width due to selective 
absorption of the continuum is quite reasonable, as long as 
two basic conditions are met: 1) there must be enough dust 
present to absorb a substantial fraction of the continuum. 2) 
this dust must be pre-dominantly located in dense neutral 
gas, so that the Lya photons are shielded from absorption. 
The latter condition is discussed in H4.5I above, and so we 
turn to discussing the first condition. 

The Lya escape fraction depends only weakly on the 
overaU dust content of the galaxy. In Fig[TSl ft^°' ~ 0.2-0.9 
over a wide range of parameters, with the escape fraction de- 
creasing with the number of surface scatters A^o and bulk gas 
motion Vc- On the other hand, because continuum photons 
are not shielded from dust by resonant scattering, they see 
the full optical depth of dust absorption, and very approxi- 
mately, /c*™ ~ exp(— Ta). A significant boost in the equiva- 
lent width therefore requires that r'' ^ 1. If the average HI 
column density across the region is (A'^) , then r'* ~ 1 requires 
(JI21 ~ i/{N2i). A damped Lya type system with an aver- 
age column density (TV) ~ 10^^ cm~^ would require a dust 
extinction cross-section per hydrogen of a'^ ~ lO"'^'^ err? /H, 
which roughly corresponds to a metalicity of ~ 1/10 so- 
lar. For sources in which these differential radiative transfer 
effects are taking place, the equivalent width should statis- 
tically have a positive correlation with the FIR flux. This 
correlation could potentially break down in more developed 
galaxies at lower redshifts, where the Lya shielding effect 
of the HI can be broken by higher gas speeds in deeper po- 
tential wells (rendering clumps optically thinner in Lya ), 
as well as the build- up/survival of dust in low density inter- 
clump gas. 



5.2 Multi-phase Outflows 

We turn to discussing the effect of a multi-phase gas out- 
flow (or inflow) on the Lya e mission line profile. The ef - 
fects of an outflowing sh e ll (e.g. iTenorio-Taele et al I il999t) ; 
lAhn. Lee, fc Led (|2003t l: lAhn l (120041) ') ar id a H ubble like 
expansion of a uniform gas sphere (e.g. iLoeb fc RvbickTI 



* However, the total dust absorption optical depth across the 
entire galaxy (many clumps) can be significantly greater than 
unity. 
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Figure 18. Lya /Continuum Escape Fraction Ratio The 

LyCK to continuum escape fraction ratio, /(f"^"//c'™> as a function 
of the total dust absorption optical depth r^ji = ^d^^— 21-^21, as- 
suming the temperature of the neutral phase is ~ lO^K. From top 
to bottom, A/q is (1, 4, 10), for bulk gas motions of 50kms~^ (solid 
lines), and 250kms~^ (dashed lines) respectively. For these pa- 
rameters, the Lya escape fractions are fc^°' = (0.94,0.68,0.38), 
and f^^" = (0.81,0.43,0.18) respectively, independent of the to- 
tal dust optical depth Ta- Most of the EW boost comes from the 
low escape fractions for continuum photons under optically thick 
conditions; very approximately, /e*™ ~ exp(— r^). 




200 



V [km/s] 



'S 






3 






i (arbitrar 


f =0.13 J \ 

esc >^^*-«w-^ \ 




3 







-1600 -1400 -1200 -1000 -800 -600 -400 -200 



V [km/s] 

Figure 19. Outflowing Shell with Gaps The normalized emis- 
sion profile as a function of the frequency shift off line center, 
for an outflowing shell with speed V-^ = 200 km/s. Two dust 
amounts were simulated: cr" = (thin line, filled in) and = 1 

(thick line). In all simulations, we used a gas temperature T4 = 1. 
For high redshift galaxies, the blue side of the profile would be 
quenched by IGM absorption. The spike of photons that escape 
at line center is easily scattered out of the line of sight by a small 
{N 1^ 5 X lO'^^ cm~^) intervening column density of HI, and thus 
not likely to be observed. 



Jl999l) : IZheng fc Miralda-Escudel (120021) ) on the Lya emis- 
sion line is a well studied problem. In both the expand- 
ing shell and expanding sphere scenarios, the generic effect 
is that a characteristic outflow speed V-^ produces a red- 
shifted emission peak at V'^"^^ ~ —V^, with an asymmetric 
shape that has a longer tail on the red side of the peak. 
The peak comprises photons that reflect off the far side of 
the expanding gas, which Doppler shifts the frequency by 
~ —V^ (see H3.4II . However, in order for these singly back- 
scattered photons to escape, the intervening gas column den- 
sity must be small enough for the photons to be transmitted, 
rather than be reflected a second time. In the rest frame of 
the near-side shell, the singly back-scattered photons have 
a frequency shift V'^'^^'' ~ -2Vf . In we showed 

that a non-negligible amount of Lya photons will be trans- 
mitted through a slab ff iVai < N^^ ^ 0.05[l/2]^[cr^2i]"^ 
Setting V ~ -2V^, we see that an outflow with speed of 
200 km/s will only allow a non-negligible amount of singly 
back-scattered photons to be transmitted if the intervening 
column density is A'^ ^ 2 x 10^''cm~'^. Observational esti- 
mates of the column density in galactic winds often exceed 
this, yet Lya is still often seen. 

The main distinguishing feature of a multiphase out- 
flow is that it allows photons of any frequency to escape 
even when the intervening gas column depth is very large, 
N ^ N^^. As in the homogeneous gas outflow models with 



smaller column densities TV < A'^''', we find that for mul- 
tiphase outflows, the emission peak is redshifted by ~ V-^ . 
However, emission is still detectible even when 2> A''^', 
as expected. 

In particular, we investigated the emission profile for 
two basic types of multiphase outflow geometries: an out- 
flowing shell with holes ( H4.2..'i^ and an outflowing ensemble 
of gas clumps modeled with the Random Surfaces geometry 
( H4.2.2t . In both cases, all surfaces were given a radially ve- 
locity with constant speed. This choice is meant to reflect 
galactic winds, where the gas reaches the asymptotic wind 
speed quickly. We placed a source of line center photons in 
the center of the region. Since the regime of optically thick 
gas has been given the least attention, we assume that the 
extreme case holds, where none of the photons penetrate 
through the gas. In this limit the surface scattering approx- 
imations of H3.3I apply in the rest frame of the scattering 
surface. The kinematics of Doppler shifts in and out of a 
moving surface can be found in H3.4I In order to distinguish 
the effects of outflow from the effects of random bulk gas mo- 
tion, we assume that there is no random bulk motion, so that 
each gas surface has an exactly radial velocity, — V^f. 
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5.2.1 Outflowing Shell with Holes 

For photons perpendicularly incident on the inner side of the 
shell, the net frequency shift induced by a moving scatter- 
ing surface is IS.V = —\V^ (see H3.4.2II . Since most photons 
escape after encountering either zero or one scattering sur- 
face, the two dominant peaks in the emission profile will be 
at V(o) = and V(i) « - fV^-''. The maximum frequency shift 
after encountering one scattering surface is AVmax = —2V^, 
and so the emission peak at V(i) should have a sharp cut-off 
at V = -2^-''. These emission peak features are verified in 
Figure [T^ The series of secondary peaks are at integer mul- 
tiples of V(i), but tend to become blended together to form 
an long red-side tail to the profile. A possible exception is 
the third peak at V(2) = 2V(i) = —^V^ (composed mostly 
of photons that encounter exactly two scattering surfaces 
before escaping), which can also be prominent in the profile. 

5.2.2 Outflowing Clumps 

For a simple model of outfiowing gas clumps we used the 
Random Surfaces geometry, which is described in detail in 
H4.2.2I In Figure 1201 we show how the profile varies with the 
covering factor, for dust free gas ct" = and for a Milky 
Way type dust content (j'L2i = 1. As for the case of an 
outflowing shell with holes, the inclusion of dust suppresses 
photons which have redshifted far from line center (so the 
HI no longer shields them from the dust), and sharpens the 
line profile. To provide some insight into how a clumpy out- 
flow causes a redshift in the emission line, in Figure ITTI we 
show how the photon radius, frequency, and cloud absorp- 
tion albedo vary as a function of n^s ■ 
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Figure 20. Outfiowing Clumps The normalized emission pro- 
file as a function of the velocity shift from line center. The thick 
lines are cr° ji = 1 while the thin lines (filled in) are dust free, 
it" = 0. The gas temperature and outflow speed are the same as 
in Figure [T^ From top to bottom, each panel is a different cov- 
ering factor: fc = 1, 3, 5. The escape fractions for the (t° ji = f 
simulations are indicated. The delta function emission spike at 
V = is composed of all the photons that escape freely without 
striking a clump. As noted in Figure IT9l exact line center pho- 
tons are likely to scattered out of the line of sight before being 
observed. 



6 CONCLUSIONS 

Our main technical, radiative transfer results are: 

• With the aid of Monte-Carlo simulations, we study the 
scattering properties of Lya photons incident on an opaque, 
dusty, moving cloud. We derive fitting formulas for the ab- 
sorption probability, frequency and angular redistribution 
functions of incident photons. 

• These formulas can be incorporated into radiative 
transfer codes, affording a vast computational speed up, and 
making feasible otherwise intractable calculations. 

• Analytically, a multiphase gas geometry can be accu- 
rately characterized by a single number, No, the number of 
surface scatters in the absence of absorption. Other factors — 
such as the cloud radii distribution for fixed A/o — are gener- 
ally unimportant. 

• We derive analytic formulas for the Lya escape fraction 
and line widths. 

• Several archetypal geometries are explored: randomly 
placed spherical clouds, randomly placed surfaces (an ab- 
straction of the prior geometry), a shell with holes, and an 
open-ended, cylindrical cavity. 

• Constant speed, radial, outflows are analyzed for bro- 
ken shells and random surfaces. The red-shifted peaks and 
widths are connected to the geometry and outflow speed. 



Our main results of direct observational relevance are: 

• Lya can escape from multi-phase dusty galaxies for HI 
column densities where it would be strongly quenched in a 
single-phase medium. 

• If most of the dust resides in a neutral phase which 
is optically thick to Lya, the Lya equivalent width can be 
strongly enhanced: while Lya photons typically scatter off 
such surfaces (which shield the dust), continuum photons 
penetrate inwards and are preferentially absorbed. 

• When the characteristic bulk gas speed exceeds ~ 
100 km/s, the Lya line width is dominated by the gas mo- 
tion, and resonant scattering frequency redistribution is sub- 
dominant effect. 

• Multiphase outflows generically produce Lya line pro- 
files that have the characteristic asymmetric shape seen in 
many starburst galaxies and Lya emitters. 

• Multiphase outflows can produce line widths several 
times larger than the actual outflow speed. 

The ISM of galaxies at both low and high redshift is 
almost certain to be both dusty and multi-phase: metal and 
dust production begin very early, given the short lifetime of 
massive stars, and thermal instability is almost inevitable 
under galactic conditions. Nonetheless, despite an extensive 
literature, to the best of our knowledge this is the first de- 



© 0000 RAS, MNRAS 000, 000-000 



24 Hansen & Oh 




Figure 21. Random Surfaces Bounce History Top: The ra- 
dial position r (in units of the moan free path between surfaces 
Ljnfp), plotted against the number of surface scatters riss- The 
thin line is a representative single photon history, the thick solid 
line is the average over an ensemble of photons. Middle: The fre- 
quency shift V (in units of the outflow speed Vf), plotted against 
riss- The thin solid line is a representative photon history, the 
thick solid line is the average history. Bottom: The Lya surface 
absorption albedo is plotted against riss, for = 200km/s 
and a gas temperature T4 = 1. We use the approximation 
e f3/^(x), which breaks down when e 1. The thin solid line 
is a representative photon history for /3 = 10~*. The thick solid 
lines are the average histories for several different values of the 
dust content (3. From lightest to darkest, /3 = 10~^°, 10"'', 10~*. 



tailed numerical study of resonance-line radiative transfer 
in a multi-phase dusty medium. The ground is surprisingly 
rich, and many future applications are envisaged! 
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APPENDIX A: TESTS OF THE CODE 

In this Appendix we show various tests of our Monte Carlo 
code. First, we test the code against known analytic solu- 
tions for optically thick slabs. Second, we compare the ac- 
celeration scheme described in >13.1.1l to exact simulations. 



Al Comparison to Analytic Solutions 

We tested the exact Monte Carlo code against known an- 
alytic solutions for very opaque slabs with a source at the 
mid-plane. First, we compared the emission frequency profile 
when the slab is pure HI, so there is no absorption. Second, 
we compared the escape fractions when the slab contains a 
small amount of absorbing dust, where the absorption cross- 
section of the dust is frequency independent. 

For an optically thick (arho > 10^), uniform slab with- 
out dust, where thq is the center-to-surfa. ce hydrogen optical 
depth at line center^ . iHarringtonl il973t) derived the mean 



^ In many pap ers on this subject, e.g. IHarringtonl ll973^ and 
iNeufeldl il990l) . "ro" refers to the mean optical depth, while 
we use the line center optical depth. In our notation, rf>{x)To = 
^{x)TfiQ where </>(x) is the normalized Voigt profile, which means 
To = \/TTTho when a < 1. 
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Figure Al. Spectra Emergent from Slab, Central Source. 

The surface mean intensity J{x) as a function of frequency, for a 
central source of Lya photons in an optically thick slab. For three 
different values of ar/iQ (labeled in the plot), the analytic surface 
intensity (grey line) is compared to exact Monte-Carlo simula- 
tions (circles). The optical depth is fixed at Tf^Q = 10®. Three dif- 
ferent gas temperatures are used: T = lOK, lO'^ K, 10'' K, which 



corresponds to ar^o = lO'* '^^, lO'^ ' 



10^ 



respectively. The 



relative error for each frequency bin is approximately l/^/Ntin, 
where Ni,i„ is the number of photons in the bin. 



intensity emission spectrum, for line center pho tons which 
are injected at the slab center (see iNeufeldl il990l) for various 
extensions) 



J{±ThO,x) 



,1/2 



2471^/2 aTho 



cosh 



J/2 



541/2 



arho 



(Al) 



where J is the mean intensity and ro is the line center hy- 
drogen optical depth from center-to-surface. Figure lXTl com- 
pares the simulation to the formula for slabs at three differ- 
ent temperatures. 

For an optically thick slab {arho > W^) with a ce nter- 
to-surface absorption optical depth Tn . INeufeldl il990l) de- 
rived the exact escape fraction for a central source, in the 
limit (aTho)^^^ S> Ta, as well as several approximation for- 
mulae. In particular, the escape fraction is well approxi- 
mated by 



/c ~ 1/ cosh 



3I/2 



{arho 



,1/3, 



(A2) 



_^5/12^ 

where is an order unity fitting parameter. INeufeldl (Il99(il 
found that the choice = 0.525 gives a good approxima- 
tion to the exact analytic escape fraction. For comparison, 
we have included in Figure I A2l the escape fractions found by 
iHummer fc Kunasd il980 ^ using numeric al integration tech- 
niques and the escape fractions found bv lAhn. Lee, fc Lee I 
l|2000|) using a Monte Carlo simulation similar to ours. We 
find that the choice of fitting parameter ( = 0.5 gives 
a good approximation to our Monte Carlo results. Both 
lAhn. Lee, fc Lee I ^00(t) and our Monte Carlo simulations 
show slightly more absorption than the analytic formula 
from INeufeldl il990l) . The analytic treatments assume the 
Lorentzian wing profile all the way down to 2; = 0, neglect- 
ing the Gaussian core. This will underestimate the number 
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Figure A2. Dusty Slab Escape Fraction, Central Source. 

The escape fraction from the middle of an optically thick, dusty 
slab is compared for several analytic and numerical methods. As 
has been shown analytically, the escape fraction depends mainly 
upon the combination (aT^o)^^'^'^a- Including dust scattering, us- 
ing e^j = 0.5 and g^j = 0.6, did not significantly affect the escape 
fraction. 



of scatters spent in the core. Although the absorption prob- 
ability per interaction is small in the core, neglecting core 
bounces will cause a slight underestimate of the overall ab- 
sorption probability. 
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Figure A3. Accelerated Monte Carlo Comparison The net 

surface absorption probability ec is shown as a function of the 
incident absorption probability per scatter e, for Lya photons 
incident on a semi-infinite plane of dusty HI at 10^ K. The dark 
symbols are simulations that use the accelerated Monte Carlo 
scheme while the white symbols are exact simulations. The circles, 
squares, and diamonds are for incident frequencies Xi = 1, 5, 20, 
respectively. The dust has an absorption albedo = 0-5 and 
a scattering asymmetry parameter g^j = 0.5. For each incident 
frequency, simulations for three values of a° were run. For each 
frequency, from left to right (or equivalently, bottom to top), the 



dust values arc cr" , 



: 0.01, 1, 100. 



A2 Testing the Acceleration Scheme 

We tested the acceleration scheme, described in H3.1.1I by 
comparing the surface absorption probability against ex- 
act simulations. Exact simulations are computationally ex- 
pensive, and so we could only test the acceleration scheme 
against a handful of exact cases. As shown in Figure lXsl the 
acceleration is quite accurate, even for initial frequencies in 
the line core. 



APPENDIX B: SURFACE Lya FREQUENCY 
REDISTRIBUTION FORMULA 

In this appendix we first show that R{xi, x; a) given in eqn. 
13911 has a unit norm over x, as claimed. Second, we outline 
the steps used to derive eqn. Mill , the generating function 
for R{xi,x; a). 

To integrate R{x,x;a) over x G (—00,00), first change 
variables to u = x^ — x^. Then the integral becomes 



f 



Ax R(xi, x; a) 



(Bl) 



The integral over the functional form (^ + it'^) ^ is a stan- 
dard integral: 

du — 7; = —7= arctan | — ^ | . (B2) 

Use of this formula in eqn. (IBH shows that the integral over 
X equals one, and so R{xi,x; a) is correctly normalized over 
X as advertised. 



To randomly generate exiting frequencies x that obey 
the probability distribution R(xi,x\ a), we use the transfor- 
mation method fr ress et alibl992.) . First, select a random 
univariate u G [0. 11 and set 



da;' R{xi, x \ a) = F{x) 



(B3) 



The frequency x is then given by functional inversion x = 
F~^{u). As above, the integral is best carried out by chang- 
ing variables to u = a;^ — if. Then F{x) is given by eqn. 
IBll l. except that the upper limit is "a;'' — if" rather than 
"00". By using equation <B2ll to complete the integration, 
we find that 



F{x) — — I arctan 

TT V 



+ 2 



(B4) 



The functional inversion gives the randomly drawn exiting 
frequency x: 

,1/3 



ii tan {tvu)) 



(B5) 



APPENDIX C: 1-D TRANSFER FOR AN 
ARBITRARY SCATTERING ASYMMETRY 
PARAMETER 

The escape fraction for arbitrary g can be approximated 
by the g = G formula, eqn. 15911 . as we demonstrate next. 
Define n* to be the average number of interactions required 
for a 50% chance of back-scattering. For interactions with 
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scattering parameter g, the probability of a forward scatter 
is (1 + g)/2 and that of a back-scatter is (1 — (?) /2. Therefore 
n* is defined to satisfy 



E 



n = 

which leads to 



1+5 



! + <? 



or equivalently 

ln(l + . 



n = 1 - 



In 2 



(CI) 



(C2) 



(C3) 



Every n* interactions acts like a single g = interaction. The 
probability of absorption after n* interactions is 1 — (1 — e)" . 
Consequently, the escape fraction is approximately the same 
as the g = G formula, eqn. I59II . with a re-scaled No and 
absorption albedo, and e*, given by 

AAo = Afo/n (C4) 
= l-(l-e)"*. (C5) 
The approximate escape fraction for arbitrary g is 
/e = l/cosh(^\/F) (C6) 
with 
Y 



= 2^(l-(l-e, 



(C7) 



where Ao ~ Noig) is calculated for the given value of g 
and n* is a function of g through eqn. IIC3II . Note that for 
1-D transfer when y = 0, then Afoig = 0) = 5(1"^ + 2r), 
which we have verified with simulations. If n* < A/o then 
there is enough scatterings for this approximation to hold. 
If en* <^ 1 also holds, then e* ~ n*e and the rescaling 
leaves the product eAo unchanged. In this limit, the escape 
fraction given in eqn. 15911 is valid for any type of scattering, 
and represents, therefore, the generic escape fraction form 
for any type of "random walk" photon transfer. On the other 
hand, if n* ^ Ao then this approximation breaks down, and 
the trajectory of the photon is more accurately characterized 
by straight line motion, with negligible back-scattering, eqn. 
15811 . As shown in Figure IHTI the approximation that the 
escape fraction is given by eqn. 15911 works well when /c ^ 
1%, and gives a decent order of magnitude estimate when 
/e is lower (such cases are observationally unimportant). 




eN„(g) 

Figure CI. 1-D Escape Fractions, Arbitrary g The sim- 
ulated escape fraction from the middle of a 1-D finite line, for 
scattering with several values of g, over a range of albedos e. 
In all simulations, the center to edge extinction (scattering + 
absorption) optical depth is constant, t = 10. Three different 
values of g were simulated (circles); from darkest to lightest 
g = 0, 1/2, —1/3. We computed Afo{g) for each value of g: 
jVoiO) = 60.0, ^"0(1/2) = 35.0, and A/'o(-l/3) = 76.3. For each 
value of g, there are five different values of e; from left to right, 
e = 0.01, 0.1, 0.3, 0.5, 0.8. The line is given by eqn. IS^ . 
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